% Prepare Output Matrices
% 1. Errors,Forecasts and RSS for Kernel Estimator
KERNEL.errors   = NaN(nfore,phi_gridsize,H_gridsize,hmax,NN_VAR);
KERNEL.forecast = NaN(nfore,phi_gridsize,H_gridsize,hmax,NN_VAR);
KERNEL.rss      = NaN(nfore,phi_gridsize,H_gridsize,NN_VAR);
% 2. Errors,Forecasts and RSS for Constant Coefficient VAR
FULL.errors   = NaN(nfore,phi_gridsize,hmax,NN_VAR);
FULL.forecast = NaN(nfore,phi_gridsize,hmax,NN_VAR);
FULL.rss      = NaN(nfore,phi_gridsize,NN_VAR);
% 3. Errors,Forecasts, RSS for RWalk forecast
RW.errors       = NaN(nfore,hmax,NN_VAR);
RW.forecast     = NaN(nfore,hmax,NN_VAR);
RW.rss          = NaN(nfore,NN_VAR);
% 4. Errors,Forecasts, RSS for Small 3 variate VAR
VAR_3.errors      = NaN(nfore,hmax,NN_VAR);
VAR_3.forecast    = NaN(nfore,hmax,NN_VAR);
VAR_3.rss         = NaN(nfore,NN_VAR);
% 5. Errors,Forecasts, RSS for N variate VAR
VAR_NN.errors      = NaN(nfore,hmax,NN_VAR);
VAR_NN.forecast    = NaN(nfore,hmax,NN_VAR);
VAR_NN.rss         = NaN(nfore,NN_VAR);
% 6. ACTUAL VALUES
ACTUALS            = NaN(nfore,hmax,NN_VAR);

for tt=esti_start:esti_end;
    if rem(tt-esti_start+1,1)==0
        disp('***************************************')
        disp([num2str(tt-esti_start+1),' out of ', num2str(esti_end-esti_start+1)]);
        disp('***************************************')
    end
    count_tt=tt-esti_start+1;
    %-- forecast window -- %
    YY_actual    = dataselect(tt+1:tt+hmax,:);
    % -- estimation window -- %
    YY_tt        = dataselect(1:tt,:);
    YY_tt_toselect = datamat(1:tt,:);		 % to be used to optimize over phi
    T_est        = size(YY_tt,1)-lags;
    constvec     = ones(T_est,1);
    % -- time t data (includes lags)-- %
    yt       = YY_tt(size(YY_tt,1)-lags+1:size(YY_tt,1),:);
    %
    % First get lags of the variables: these are the same for all the
    % models so it is more efficient to compute them outside of the routines
    [y_l,x_l] = matrix_var(YY_tt,lags,T_est+lags,m);
    % Fit single equation models to each equation to get an estimate of the covariance matrix
    % This also remains the same across models
    sigu_ar = zeros(1,m);
    for ll=1:m;
        y_l_ar = y_l(:,ll);
        pos_ar = seqa(ll,m,lags);
        x_l_ar = [x_l(:,pos_ar') constvec];
        Proj_x = x_l_ar*((x_l_ar'*x_l_ar)\x_l_ar')*y_l_ar;
        sigu_ar(1,ll)=(y_l_ar-Proj_x)'*(y_l_ar-Proj_x)/T_est;
    end;
    
    poslags = [vec([1:m:lags*m; 2:m:lags*m; 3:m:lags*m])' m*lags+1];
    %% compute matrices for VAR estimation that are constant for all models
    % 1. matrix Omega
    sigmas=sigu_ar.^.5;
    Omega=[kron(diag(1:lags),diag(sigmas)) zeros(m*lags,1);        zeros(1,m*lags+1)];
    Omega(end,end)=1/(1e+20);%//very lose prior on the constant like banbura et al.
    % The matrix OmegaprimeOmega is the "Ridge" term that arises when implementing the dummy
    % observation approach in Banbura et al.
    OmegaprimeOmega = Omega'*Omega;
    % save in cells the VAR coefficients for TVP-VAR (one for each point in the Phi/H grid)
    slopes           = cell(length(phi_grid),length(H_grid));
    % save in cells the VAR coefficients for ConstantParams-VAR (one for each point in the Phi)
    slopes_nokernel  = cell(1,length(H_grid));
    tic
    d1               = [diag(gamma_mean);zeros(m*(lags-1)+1,m)];
    close all;
    conta            = 0;

    for posgrid_phi=1:phi_gridsize
        big_X=x_l'*x_l;
        big_Y=x_l'*y_l;
        phi     = phi_grid(posgrid_phi);
    % Coefficients for large VAR with Constant Parameters       
        b_full  = (big_X+OmegaprimeOmega/(phi^2))\(big_Y+(OmegaprimeOmega/(phi^2))*d1);       
        slopes_nokernel{posgrid_phi,1}=b_full;
        [y_full,rss_full]                                = var_prior(posgrid_phi,1,lags,hmax,m,T_est,y_l,x_l,slopes_nokernel);
        %%
        % 1. FULL SAMPLE VAR WITH PRIORS
        FULL.errors(count_tt,posgrid_phi,:,:)  = y_full(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
        FULL.forecast(count_tt,posgrid_phi,:,:)= y_full(:,1:NN_VAR);
        FULL.rss(count_tt,posgrid_phi,:)       = rss_full(1:NN_VAR);
        for posgrid_H=1:H_gridsize
            % compute the kernel weights
            if strcmp('rect', kernshape)
                Wt   = w_const(T_est,H_grid(posgrid_H));
            elseif strcmp('gauss', kernshape)
                H = T_est.^H_grid(posgrid_H);
                % when using gauss , third input must be 1 
                Wt = w_gauss(T_est,T_est,H,1);
            elseif strcmp('ewma', kernshape)
                H = T_est.^H_grid(posgrid_H);
                Wt = w_ewma(T_est,H_grid(posgrid_H));
            end;
% HERE COMPUTE KERNEL BASED VAR, first compute Weight Matrix Wt, then
% perform weighted least squared estimation.
% To speed up computation approximate the number of observations
% Throw away weights lower than 1/10000
            T_approx = sum((diag(Wt))<(1/10^4));
            y_l_appr=y_l(T_approx+1:end,:);
            x_l_appr=x_l(T_approx+1:end,:);
            T_est_appr = T_est-T_approx;
            temp = diag(Wt);
            Wt_appr = diag(temp(T_approx+1:end));
            big_X_appr=x_l_appr'*Wt_appr*x_l_appr;
            big_Y_appr=x_l_appr'*Wt_appr*y_l_appr;
            phi = phi_grid(posgrid_phi);
            b  = (big_X_appr+OmegaprimeOmega/(phi^2))\(big_Y_appr+(OmegaprimeOmega/(phi^2))*d1);
            slopes{posgrid_phi,posgrid_H}=b;
            [y_kern,rss_kern]   = var_prior(posgrid_phi,posgrid_H,lags,hmax,m,T_est,y_l,x_l,slopes);
            % 2. KERNEL BASED VAR WITH PRIORS
            KERNEL.errors(count_tt,posgrid_phi,posgrid_H,:,:)  = y_kern(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
            KERNEL.forecast(count_tt,posgrid_phi,posgrid_H,:,:)= y_kern(:,1:NN_VAR);
            KERNEL.rss(count_tt,posgrid_phi,posgrid_H,:)       = rss_kern(1:NN_VAR);
        end
    end

    %% now estimate the same model with
    % 3. TRIVARIATE VAR (OLS)
    poslags                         = [vec([1:m:lags*m; 2:m:lags*m; 3:m:lags*m])' m*lags+1];
    [y_ff_3,rss_3]                  = var_ols(lags,hmax,3,T_est,y_l(:,1:3),x_l(:,poslags));
    VAR_3.errors(count_tt,:,:)      = y_ff_3(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
    VAR_3.forecast(count_tt,:,:)    = y_ff_3(:,1:NN_VAR);
    VAR_3.rss(count_tt,:)           = rss_3(1:NN_VAR);
    % 4. N VARIATE VAR (OLS)
    [y_ff_NN,rss_NN]                = var_ols(lags,hmax,m,T_est,y_l,x_l);
    VAR_NN.errors(count_tt,:,:)     = y_ff_NN(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
    VAR_NN.forecast(count_tt,:,:)   = y_ff_NN(:,1:NN_VAR);
    VAR_NN.rss(count_tt,:)          = rss_NN(1:NN_VAR);
    % 5. RANDOM WALK WITH DRIFT (OLS)
    [y_ff_RW,rss_RW]                = var_RW(hmax,m,T_est,y_l,x_l);
    RW.errors(count_tt,:,:)         = y_ff_RW(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
    RW.forecast(count_tt,:,:)       = y_ff_RW(:,1:NN_VAR);
    RW.rss(count_tt,:)              = rss_RW(1:NN_VAR);
    ACTUALS(count_tt,:,:)           = YY_actual(:,1:NN_VAR);
end

%% SAVE RESULTS
 save(FileOutput);

% nn=400;
% aa=squeeze(KERNEL.rss(nn,:,end,:));
% bb=squeeze(FULL.rss(nn,:,:));
% cc=repmat(VAR_3.rss(nn,:),phi_gridsize,1);
% dd=repmat(VAR_NN.rss(nn,:),phi_gridsize,1);
% ee=repmat(RW.rss(nn,:),phi_gridsize,1);
% lll=0;
% for ll=1:3; figure(ll); lll=lll+[aa(:,ll) bb(:,ll) cc(:,ll) dd(:,ll) ee(:,ll)];plot([aa(:,ll) bb(:,ll) cc(:,ll) dd(:,ll) ee(:,ll)]); legend('kern','full','var3','VARNN','RW');end
% close all; plot(lll)