function [OUT] = VAR_estimate_kernel_6(data_l,m,d1,kernshape,esti_start,esti_end,phi_grid,H_grid,phi_gridsize,H_gridsize,hmax,NN_VAR,lags,const)
% uguale a VAR_estimate_kernel_4 ma stima OmegaprimeOmega a ogni passo

d1(2) = 0;
% do no standardize data
data    = data_l; 
std_mat = 1;
med_mat = 0;  
nfore        = esti_end-esti_start+1;
ncoef        = lags*m^2+m*const;
% ALLOCATE DIMENSION FOR OUTPUT
for shapes=1:length(kernshape)
    OUT.(kernshape{shapes}).errors_sel    = NaN(nfore,hmax,NN_VAR);
    OUT.(kernshape{shapes}).errors_ave1   = NaN(nfore,hmax,NN_VAR);
    OUT.(kernshape{shapes}).errors_ave2   = NaN(nfore,hmax,NN_VAR);
    OUT.(kernshape{shapes}).forecast      = NaN(nfore,hmax,NN_VAR);
    OUT.(kernshape{shapes}).coeffs        = NaN(nfore,ncoef);
end
clear shapes;

for tt=esti_start:esti_end;
%     if rem(tt,20)==0
%         disp(['Estimating kernel based VAR t=',num2str(tt)]);
%     end
    count_tt=tt-esti_start+1;
%    YY_actual    = data(tt+1:tt+hmax,:);
    YY_actual    = data_l(tt+1:tt+hmax,:);
    
    [y_l,x_l_temp]  = matrix_var(data(1:tt,:),lags,tt,m);
    [~,~,OmegaprimeOmegatemp] = get_Omega(tt-lags,m,lags,y_l,x_l_temp);
    OmegaprimeOmega = OmegaprimeOmegatemp(1:end-1+const,1:end-1+const);
    x_l   = x_l_temp(:,1:end-1+const);
    T_est = size(x_l,1);
    for shapes=1:length(kernshape)
        for posgrid_H=1:H_gridsize
            if strcmp('rect', kernshape{shapes})
                Wt   = w_const(T_est,H_grid(posgrid_H));
            elseif strcmp('gauss', kernshape{shapes})
                H = T_est.^H_grid(posgrid_H);
                Wt = w_gauss(T_est,T_est,H,1);
            elseif strcmp('ewma', kernshape{shapes})
                Wt = w_ewma(T_est,H_grid(posgrid_H));
            end;
            for posgrid_phi=1:phi_gridsize
                % approximate the number of observations to further speed up
                % computation: throw away weights lower than 1/10000
                T_approx = 0;% sum((diag(Wt))<(1/10^4));
                y_l_appr = y_l(T_approx+1:end,:);
                x_l_appr = x_l(T_approx+1:end,:);
                temp     = diag(Wt);
                Wt_appr  = sparse(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);
                TT = T_est;
                %
                %                 %      now perform recursive forecasts
                y_f = y_l(TT,:);
                x_f = x_l(TT,:);
                y_ff = NaN(hmax,m);
                for llags = 1:hmax
                    x_f  = [y_f x_f(1,1:m*lags-m) 1];
                    xff  = x_f(:,1:end-(1-const));
                    y_f  = xff*b;
                    y_ff(llags,:) = y_f.*std_mat(1,:)+med_mat(1,:);
                end;
                %                 % COLLECT OUTPUT
                TEMPOUT.errors.(kernshape{shapes})(count_tt,posgrid_phi,posgrid_H,:,:)  = y_ff(:,1:NN_VAR)-YY_actual(:,1:NN_VAR);
                if sum(isnan(y_ff(:,1:NN_VAR)))==0
                else
                    keyboard
                end
                    
                TEMPOUT.forecast.(kernshape{shapes})(count_tt,posgrid_phi,posgrid_H,:,:)= y_ff(:,1:NN_VAR);
                TEMPOUT.coeffs.(kernshape{shapes})(count_tt,posgrid_phi,posgrid_H,:)    = b(:);
            end
        end
    end
end

%% Compute model weights.
% %  Start from computing squared errors
disp('compute model weights');
[~,stad,~] = normaliz(data_l);
window = 40;
for shapes=1:length(kernshape)
    model_space = phi_gridsize*H_gridsize;
    ERR  = reshape(TEMPOUT.errors.(kernshape{shapes}),nfore,model_space,hmax,NN_VAR);
    FORE = reshape(TEMPOUT.forecast.(kernshape{shapes}),nfore,model_space,hmax,NN_VAR);
    % COEF = reshape(TEMPOUT.coeffs.(kernshape{shapes}),nfore,model_space,ncoef);
    % start from positing equal weights
    p_mod_stand       = ones(nfore,model_space,hmax)./(model_space);
    p_mod_temp        = p_mod_stand;

    for hh=1:hmax
        disp(['horizon=',num2str(hh)])
        for tt=hh+1:nfore;
            for pos=1:model_space
    % you need to "normalize" squared errors, since you are averaging
    % across the cross-section    
            if tt<=window+hmax
                if tt-hh==1
                    err = (squeeze(ERR(1:tt-hh,pos,hh,:)))';     % only see the error h periods before
                else
                    err = (squeeze(ERR(1:tt-hh,pos,hh,:)));      % only see the error h periods before
                end
            else
                err = (squeeze(ERR(tt-window-hh:tt-hh,pos,hh,:)));      % only see the error h periods before
            end
                std_sqerr  = bsxfun(@rdivide,err.^2,stad(1,1:NN_VAR).^2);
%               mean_sqerr = mean(std_sqerr,2); %this is the mean in the cross section, so this is a (t-hh)x1 vector
                mean_sqerr = std_sqerr; %this is the mean in the cross section, so this is a (t-hh)x1 vector
                p_mod_temp(tt,pos,hh)  = mean(1./((mean(mean_sqerr)).^.5),2);
                %            p_mod_ewma_temp(tt,hh,pos)  = 1./(mean(mat_err_sq)).^.5;
            end
            %    normalize weights
            p_mod_stand(tt,:,hh)       = p_mod_temp(tt,:,hh)./(sum(p_mod_temp(tt,:,hh),2));
        end
    end
end
    

    % now take EWMA exponentials of p_mod
    %     p_mod_disc = NaN(tt,model_space);
    %     p_mod_disc(1,:)=log(p_mod_stand(1,:));% initial condition
    %
    %     for tt=2:nfore
    %         vec_sum = ones(1,tt);
    %         W_t = w_ewma(tt,.99);
    %         p_mod_disc(tt,:)=vec_sum*(W_t*log(p_mod_stand(1:tt,:))); %questo fa la somma dei log square fore errors, discounted con EWMA
    %     end
    %
    %     p_mod_exp    =  exp(p_mod_disc);
    %     denominator  =  repmat(sum(p_mod_exp,2),1,model_space);
    %
    %     p_mod_stand         = p_mod_exp./denominator;
    for gg=1:hmax
        [~,ind_DMS(:,gg)]         = max(p_mod_stand(:,:,gg),[],2);
    end
    %     % err_opt             = squeeze(ERR(:,ind_DMS,:,:));
    % %    err_opt       =  NaN(tt,NN_VAR);
    equal_w              = p_mod_stand*0+1/model_space;
    OUT.(kernshape{shapes}).errors_ave1   =   squeeze(sum(bsxfun(@times,p_mod_stand,ERR),2));
    OUT.(kernshape{shapes}).errors_ave2   =   squeeze(sum(bsxfun(@times,equal_w,ERR),2));
    
    for ll=1:nfore
        for hh=1:hmax
            temp = ERR(ll,ind_DMS(ll,hh),hh,:);
            if sum(isnan(temp)==0)
            else
                keyboard
            end
            OUT.(kernshape{shapes}).errors_sel(ll,hh,:)   =  temp;
            %        OUT.(kernshape{shapes}).forecast(ll,hh,:) =  FORE(ll,ind_DMS(ll,hh),1:hmax,:);
            %        OUT.(kernshape{shapes}).coeffs(ll,:,:)   =  COEF(ll,ind_DMS(ll),:);
            clear temp
        end
    end
end

