% Performs the adaptive elastic net with standardization
% for a vector of hyper parameters lambda, and hyperparameter alpha = a
% uses regularization Weights.
%
% This function is needed to ensure that the standardization is done
% consistently over the regression.
% Removes the intercept before the regression and adds it after.
% In particular, it is needed for a VAR specification.

function [estimator] = adaptiveElasticNet(X,Z,lambda,a,Weights)
% find sizes first
if size(X,1)==size(Z,1)-1
    intercept = 1;
elseif size(X,1)==size(Z,1)
    intercept = 0;
else
    intercept = 1;
    factors = size(Z,1) - size(X,1) -1;
end

% standardarize independent data to ensure that no variable is represented
% too strong
Xa = (X-mean(X,2))./std(X')'; 
Za = Z;
% standardize dependent variable without intercept.
Za(1+intercept:end,:) = (Z(1+intercept:end,:)-mean(Z(1+intercept:end,:),2))./std(Z(1+intercept:end,:)')';

% if no weights, use regularization on all but the factors and intercept
if nargin > 4
    weights = Weights;
else
    weightTemp = ones(size(X,1),size(Z,1));
    if intercept >= 1
        weightTemp(:,1:intercept) = inf; %Factors
    end
    weights = (weightTemp(:)).^(-1);
end

% exclude variables with weight infinity to speed up the optimization
auxVec = 1:length(weights);
exclude = auxVec(weights == inf)';
weights(weights == inf) = 0;

% define options and run glmnet
opts = struct('alpha', a,'intr',0,'lambda',lambda,'standardize',0,'penalty_factor',weights,'exclude',exclude);
fit = glmnet(sparse( kron(Za',speye(size(X,1)))), Xa(:), [], opts);

% refit the estimator towards the non-standardized variables
for ii=1:length(lambda)
    estimator(:,:,ii) = reshape( fit.beta(:,ii), size(X,1), size(Z,1)) ;  % add weight if needed: (  1+(1-a)*lambda(ii) )*
    estimator(:,1+intercept:end,ii) = diag(std(X'))*estimator(:,1+intercept:end,ii)/diag(std(Z(1+intercept:end,:)'));
    estimator(:,1:intercept,ii) = diag(std(X'))*estimator(:,1:intercept,ii) + mean(X,2) -estimator(:,1+intercept:end,ii)*mean(Z(1+intercept:end,:),2);
end
estimator = squeeze(estimator);
end

