% This is an executable file which runs the various types of forecasting
% models (simple models, model averages, factor model, var model) and 
% stores them in matlab data structures PITSgdp1, PITSgdp4, PITSpgdp1, PITSpgdp4

clc; clear all; close all;

param.R = 40;
param.pmax = 4;
param.nkeep = 10000;
qn = floor(param.R^(1/4));
h = [1 4];

for hzn = h
    % Put the data in appropriate objects
    data = getData(hzn);

    param.tgdp = length(data.gdp);
    param.tpgdp = length(data.pgdp);
    param.tvar = length(data.pgdp);
    
    dategdp = data.tseriesgdp(param.R+1:end,:);
    datepgdp = data.tseriespgdp(param.R+1:end,:);

    vnamegdp = data.vnamegdp;
    vnamepgdp = data.vnamepgdp;
    
    ygdp = data.gdp;
    ypgdp = data.pgdp;
    xgdp = data.xgdp;
    xpgdp = data.xpgdp;
    regdatagdp = data.regdatagdp;
    regdatapgdp = data.regdatapgdp;
    
    vary = data.vary;
    varx = data.varx;
    
    % Creates matrices to store the forecasts
    densforecastgdp = NaN(size(regdatagdp,2)+1,param.tgdp-param.R);
    densforecastpgdp = NaN(size(regdatapgdp,2)+1,param.tpgdp-param.R);
    scoredensforecastgdp = NaN(size(regdatagdp,2)+1,param.tgdp-param.R);
    scoredensforecastpgdp = NaN(size(regdatapgdp,2)+1,param.tpgdp-param.R);

    bayesdensforecastgdp = NaN(size(regdatagdp,2)+1,param.tgdp-param.R);
    bayesdensforecastpgdp = NaN(size(regdatapgdp,2)+1,param.tpgdp-param.R);
    scorebayesdensforecastgdp = NaN(size(regdatagdp,2)+1,param.tgdp-param.R);
    scorebayesdensforecastpgdp = NaN(size(regdatapgdp,2)+1,param.tpgdp-param.R);
    bayesmlgdp = NaN(size(regdatagdp,2)+1,param.tgdp-param.R);
    bayesmlpgdp = NaN(size(regdatapgdp,2)+1,param.tpgdp-param.R);
    
    pooldensforecastgdp = NaN(4,param.tgdp-param.R);
    pooldensforecastpgdp = NaN(4,param.tpgdp-param.R);
    scorepooldensforecastgdp = NaN(4,param.tgdp-param.R);
    scorepooldensforecastpgdp = NaN(4,param.tpgdp-param.R);

    bmamlgdp = NaN(param.nkeep,param.tgdp-param.R);
    bmamlpgdp = NaN(param.nkeep,param.tpgdp-param.R);
    bmadensforecastgdp = NaN(param.nkeep,param.tgdp-param.R);
    bmadensforecastpgdp = NaN(param.nkeep,param.tpgdp-param.R);
    scorebmadensforecastgdp = NaN(param.nkeep,param.tgdp-param.R);
    scorebmadensforecastpgdp = NaN(param.nkeep,param.tpgdp-param.R);
    
    vardensforecastgdp = NaN(1,param.tvar-param.R);
    vardensforecastpgdp = NaN(1,param.tvar-param.R);
    scorevardensforecastgdp = NaN(1,param.tvar-param.R);
    scorevardensforecastpgdp = NaN(1,param.tvar-param.R);

    factordensforecastgdp = NaN(1,param.tgdp-param.R);
    factordensforecastpgdp = NaN(1,param.tpgdp-param.R);
    scorefactordensforecastgdp = NaN(1,param.tgdp-param.R);
    scorefactordensforecastpgdp = NaN(1,param.tpgdp-param.R);
    
    % Initialize the values necessery to calibrate the g-prior for BMA
    [betapriorgdp betapriorpgdp] = prior(hzn);
    phi = 1;

    % OUTPUT GROWTH FORECASTS
    % one observation is lost to the out-of-sample forecast
    for i = param.R:(param.tgdp-1)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % AR model 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        roll_ygdp = ygdp(i-param.R+1:i,:);
        roll_xgdp = xgdp(i-param.R+1:i,:);
        
        param.l = mbic(roll_ygdp,ones(param.R,1),roll_xgdp, param.pmax);

        rgrgdp = ones(param.R-param.l+1,1);

        for j = 1:param.l
            rgrgdp = [rgrgdp roll_xgdp(param.l-j+1:end-j+1,:)];
        end
        
        [beta betaint e] = regress(roll_ygdp(param.l:end,:),rgrgdp);
        resultsgdp.beta = beta;
        resultsgdp.sige = nw(e,qn);
        
        oosgdp = [1 xgdp(i+1,:) rgrgdp(end,2:end-1)];
        densforecastgdp(1,i-param.R+1) = normcdf(ygdp(i+1,1),oosgdp*resultsgdp.beta,sqrt(resultsgdp.sige));
        scoredensforecastgdp(1,i-param.R+1) = normpdf(ygdp(i+1,1),oosgdp*resultsgdp.beta,sqrt(resultsgdp.sige));

        % Bayesian pooled models
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % BMA
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [bayesgdp.beta bayesgdp.sige] = bmareg(roll_ygdp(param.l:end,:),rgrgdp,[betapriorgdp; zeros(param.l-1,1)],phi);
        bayesdensforecastgdp(1,i-param.R+1) = normcdf(ygdp(i+1,1),oosgdp*bayesgdp.beta,sqrt(bayesgdp.sige));
        scorebayesdensforecastgdp(1,i-param.R+1) = normpdf(ygdp(i+1,1),oosgdp*bayesgdp.beta,sqrt(bayesgdp.sige));
        bayesmlgdp(1,i-param.R+1) = exp(modellikel(roll_ygdp(param.l:end,:),rgrgdp,[betapriorgdp; zeros(param.l-1,1)],phi));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % BMA via MC3
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [bmamlgdp(:,i-param.R+1), bmadensforecastgdp(:,i-param.R+1)  scorebmadensforecastgdp(:,i-param.R+1)] = covselection(ygdp(i+1,1),roll_ygdp(param.l:end,:),rgrgdp,regdatagdp(i-param.R+param.l:i,:), [betapriorgdp; zeros(param.l-1,1)],1,2,param.nkeep,[oosgdp regdatagdp(i+1,:)]);
        bmamlgdp(:,i-param.R+1) = bmamlgdp(:,i-param.R+1)./sum(bmamlgdp(:,i-param.R+1));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % factor model
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        gdpic2 = ICP(cleanNaN(regdatagdp(i-param.R+param.l:i,:)),10,8,2);
        [gdpehat,gdpFhat,gdplamhat,gdpve2] = pc(standard(cleanNaN(regdatagdp(i-param.R+param.l:i,:))),gdpic2); % get rhat factors using L'L/N normalization
        [beta betaint e] = regress(roll_ygdp(param.l:end,:),[rgrgdp gdpFhat]);
        resultsfactorgdp.beta = beta;
        resultsfactorgdp.sige = nw(e,qn);
        
        [gdpehattemp,gdpFhattemp,gdplamhattemp,gdpve2temp] = pc(standard(cleanNaN(regdatagdp(i-param.R+param.l+1:i+1,:))),gdpic2); % get rhat factors using L'L/N normalization
        oosfactorgdp = [oosgdp gdpFhattemp(end,:)];
        
        factordensforecastgdp(1,i-param.R+1) = normcdf(ygdp(i+1,1),oosfactorgdp*resultsfactorgdp.beta,sqrt(resultsfactorgdp.sige));
        scorefactordensforecastgdp(1,i-param.R+1) = normpdf(ygdp(i+1,1),oosfactorgdp*resultsfactorgdp.beta,sqrt(resultsfactorgdp.sige));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Models with Additional Regressors
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for k = 1:size(regdatagdp,2)
            roll_reggdp = regdatagdp(i-param.R+1:i,k);

            if (isfinite(regdatagdp(i-param.R+1:i+1,k)) == 1) 
                param.rl = mbic(roll_ygdp(param.l:end),rgrgdp,roll_reggdp(param.l:end), param.pmax);
                param.ml = max(param.rl, param.l);

                rgrgdpreg = ones(param.R-param.ml+1,1);

                for j = 1:param.l
                    rgrgdpreg = [rgrgdpreg roll_xgdp(param.ml-j+1:end-j+1,:)];
                end
                for j = 1:param.rl
                    rgrgdpreg = [rgrgdpreg roll_reggdp(param.ml-j+1:end-j+1,:)];
                end

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % ADL model 
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                [beta betaint e] = regress(roll_ygdp(param.ml:end,:),rgrgdpreg);
                resultsgdp.beta = beta;
                resultsgdp.sige = nw(e,qn);
                
                oosgdp = [1 xgdp(i+1,:) rgrgdpreg(end,2:2+param.l-2)  regdatagdp(i+1,k)  rgrgdpreg(end,1+param.l+1:1+param.l+1+param.rl-2)]; %-2 are there since were are trying to truncate one observation for oos forecasting.
                densforecastgdp(k+1,i-param.R+1) = normcdf(ygdp(i+1,1),oosgdp*resultsgdp.beta,sqrt(resultsgdp.sige));
                scoredensforecastgdp(k+1,i-param.R+1) = normpdf(ygdp(i+1,1),oosgdp*resultsgdp.beta,sqrt(resultsgdp.sige));

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % Bayesian pooled models
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                bhatpriorgdp = [betapriorgdp; zeros(param.l-1,1); zeros(param.rl,1)];
                [bayesgdp.beta bayesgdp.sige] = bmareg(roll_ygdp(param.ml:end,:),rgrgdpreg,bhatpriorgdp,phi);
                bayesdensforecastgdp(k+1,i-param.R+1) = normcdf(ygdp(i+1,1),oosgdp*bayesgdp.beta,sqrt(bayesgdp.sige));
                scorebayesdensforecastgdp(k+1,i-param.R+1) = normpdf(ygdp(i+1,1),oosgdp*bayesgdp.beta,sqrt(bayesgdp.sige));
                bayesmlgdp(k+1,i-param.R+1) = exp(modellikel(roll_ygdp(param.ml:end,:),rgrgdpreg,bhatpriorgdp,phi));
            end
        end
    end
    
    pooldensforecastgdp(1,:) = (sum(cleanNaN(densforecastgdp(2:end,:)'),2)./size(cleanNaN(densforecastgdp(2:end,:)'),2))';
    pooldensforecastgdp(2,:) = (sum(cleanNaN(densforecastgdp(2:end,:)').*cleanNaN(bayesmlgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlgdp(2:end,:)'),2)),2))';
    pooldensforecastgdp(3,:) = (sum(cleanNaN(bayesdensforecastgdp(2:end,:)').*cleanNaN(bayesmlgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlgdp(2:end,:)'),2)),2))';
    pooldensforecastgdp(4,:) = sum(bmamlgdp.*bmadensforecastgdp,1);
    
    scorepooldensforecastgdp(1,:) = (sum(cleanNaN(scoredensforecastgdp(2:end,:)'),2)./size(cleanNaN(scoredensforecastgdp(2:end,:)'),2))';
    scorepooldensforecastgdp(2,:) = (sum(cleanNaN(scoredensforecastgdp(2:end,:)').*cleanNaN(bayesmlgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlgdp(2:end,:)'),2)),2))';
    scorepooldensforecastgdp(3,:) = (sum(cleanNaN(scorebayesdensforecastgdp(2:end,:)').*cleanNaN(bayesmlgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlgdp(2:end,:)'),2)),2))';
    scorepooldensforecastgdp(4,:) = sum(bmamlgdp.*scorebmadensforecastgdp,1);


    % INFLATION FORECASTS
    % one observation is lost to the out-of-sample forecast
    for i = param.R:(param.tpgdp-1)
        % BVAR model
        roll_vary = data.vary(i-param.R+1:i,:);
        roll_varx = data.varx(i-param.R+1:i,:);
       
        param.varl = mbic(cleanNaN(roll_vary),ones(param.R,1),cleanNaN(roll_varx), param.pmax);

        vx = ones(param.R-param.varl+1,1);
        for j = 1:param.varl
            vx = [vx roll_varx(param.varl-j+1:end-j+1,:)];
        end
        
        oosvar = [1 varx(i+1,:) vx(end,2:end-size(varx,2))];
        vartemp = varestimate(roll_vary(param.varl:end,:),vx,vary(i+1,:),oosvar);
        
        vardensforecastgdp(1,i-param.R+1) = vartemp.pitgdp;
        vardensforecastpgdp(1,i-param.R+1) = vartemp.pitpgdp;
        scorevardensforecastgdp(1,i-param.R+1) = vartemp.scoregdp;
        scorevardensforecastpgdp(1,i-param.R+1) = vartemp.scorepgdp;
 
        % AR model
        roll_ypgdp = ypgdp(i-param.R+1:i,:);
        roll_xpgdp = xpgdp(i-param.R+1:i,:);

        param.pl = mbic(roll_ypgdp,ones(param.R,1),roll_xpgdp, param.pmax);

        rgrpgdp = ones(param.R-param.pl+1,1);

        for j = 1:param.pl
            rgrpgdp = [rgrpgdp roll_xpgdp(param.pl-j+1:end-j+1,:)];
        end

        [beta betaint e] = regress(roll_ypgdp(param.pl:end,:),rgrpgdp);
        resultspgdp.beta = beta;
        resultspgdp.sige = nw(e,qn);
  
        oospgdp = [1 xpgdp(i+1,:) rgrpgdp(end,2:end-1)];
        densforecastpgdp(1,i-param.R+1) = normcdf(ypgdp(i+1,1),oospgdp*resultspgdp.beta,sqrt(resultspgdp.sige));
        scoredensforecastpgdp(1,i-param.R+1) = normpdf(ypgdp(i+1,1),oospgdp*resultspgdp.beta,sqrt(resultspgdp.sige));
    
        % Bayesian pooled models
        % BMA
        [bayespgdp.beta bayespgdp.sige] = bmareg(roll_ypgdp(param.pl:end,:),rgrpgdp,[betapriorpgdp; zeros(param.pl-1,1)],phi);
        bayesdensforecastpgdp(1,i-param.R+1) = normcdf(ypgdp(i+1,1),oospgdp*bayespgdp.beta,sqrt(bayespgdp.sige));
        scorebayesdensforecastpgdp(1,i-param.R+1) = normpdf(ypgdp(i+1,1),oospgdp*bayespgdp.beta,sqrt(bayespgdp.sige));
        bayesmlpgdp(1,i-param.R+1) = exp(modellikel(roll_ypgdp(param.pl:end,:),rgrpgdp,[betapriorpgdp; zeros(param.pl-1,1)],phi));
        
        % BMA via MC3
        [bmamlpgdp(:,i-param.R+1), bmadensforecastpgdp(:,i-param.R+1), scorebmadensforecastpgdp(:,i-param.R+1)] = covselection(ypgdp(i+1,1), roll_ypgdp(param.pl:end,:),rgrpgdp,regdatapgdp(i-param.R+param.pl:i,:), [betapriorpgdp; zeros(param.pl-1,1)],1,2,param.nkeep,[oospgdp regdatapgdp(i+1,:)]);
        bmamlpgdp(:,i-param.R+1) = bmamlpgdp(:,i-param.R+1)./sum(bmamlpgdp(:,i-param.R+1));

        % Factor models
        pgdpic2 = ICP(cleanNaN(regdatapgdp(i-param.R+param.pl:i,:)),10,8,2);
        [pgdpehat,pgdpFhat,pgdplamhat,pgdpve2] = pc(standard(cleanNaN(regdatapgdp(i-param.R+param.pl:i,:))),pgdpic2); % get rhat factors using L'L/N normalization
        [beta betaint e] = regress(roll_ypgdp(param.pl:end,:),[rgrpgdp pgdpFhat]);
        resultsfactorpgdp.beta = beta;
        resultsfactorpgdp.sige = nw(e,qn);
        
        [pgdpehattemp,pgdpFhattemp,pgdplamhattemp,pgdpve2temp] = pc(standard(cleanNaN(regdatapgdp(i-param.R+param.pl+1:i+1,:))),pgdpic2); % get rhat factors using L'L/N normalization
        oosfactorpgdp = [oospgdp pgdpFhattemp(end,:)];
       
        factordensforecastpgdp(1,i-param.R+1) = normcdf(ypgdp(i+1,1),oosfactorpgdp*resultsfactorpgdp.beta,sqrt(resultsfactorpgdp.sige));     
        scorefactordensforecastpgdp(1,i-param.R+1) = normpdf(ypgdp(i+1,1),oosfactorpgdp*resultsfactorpgdp.beta,sqrt(resultsfactorpgdp.sige));     

        for k = 1:size(regdatapgdp,2)
            roll_regpgdp = regdatapgdp(i-param.R+1:i,k);

            if (isfinite(regdatapgdp(i-param.R+1:i+1,k)) == 1) 
                param.prl = mbic(roll_ypgdp(param.pl:end),rgrpgdp,roll_regpgdp(param.pl:end), param.pmax);
                param.pml = max(param.prl, param.pl);

                rgrpgdpreg = ones(param.R-param.pml+1,1);

                for j = 1:param.pl
                    rgrpgdpreg = [rgrpgdpreg roll_xpgdp(param.pml-j+1:end-j+1,:)];
                end
                for j = 1:param.prl
                    rgrpgdpreg = [rgrpgdpreg roll_regpgdp(param.pml-j+1:end-j+1,:)];
                end

                % ADL models
                [beta betaint e] = regress(roll_ypgdp(param.pml:end,:),rgrpgdpreg);
                resultspgdp.beta = beta;
                resultspgdp.sige = nw(e,qn);

                oospgdp = [1 xpgdp(i+1,:) rgrpgdpreg(end,2:2+param.pl-2)  regdatapgdp(i+1,k)  rgrpgdpreg(end,1+param.pl+1:1+param.pl+1+param.prl-2)]; %-2 are there since were are trying to truncate one observation for oos forecasting.
                densforecastpgdp(k+1,i-param.R+1) = normcdf(ypgdp(i+1,1),oospgdp*resultspgdp.beta,sqrt(resultspgdp.sige));
                scoredensforecastpgdp(k+1,i-param.R+1) = normpdf(ypgdp(i+1,1),oospgdp*resultspgdp.beta,sqrt(resultspgdp.sige));
                
                % Bayesian pooled models
                % BMA
                bhatpriorpgdp = [betapriorpgdp; zeros(param.pl-1,1); zeros(param.prl,1)];
                [bayespgdp.beta bayespgdp.sige] = bmareg(roll_ypgdp(param.pml:end,:),rgrpgdpreg,bhatpriorpgdp,phi);
                bayesdensforecastpgdp(k+1,i-param.R+1) = normcdf(ypgdp(i+1,1),oospgdp*bayespgdp.beta,sqrt(bayespgdp.sige));
                scorebayesdensforecastpgdp(k+1,i-param.R+1) = normpdf(ypgdp(i+1,1),oospgdp*bayespgdp.beta,sqrt(bayespgdp.sige));
                bayesmlpgdp(k+1,i-param.R+1) = exp(modellikel(roll_ypgdp(param.pml:end,:),rgrpgdpreg,bhatpriorpgdp,phi));
            end
        end
    end
    
   pooldensforecastpgdp(1,:) = (sum(cleanNaN(densforecastpgdp(2:end,:)'),2)./size(cleanNaN(densforecastpgdp(2:k+1,:)'),2))';
   pooldensforecastpgdp(2,:) = (sum(cleanNaN(densforecastpgdp(2:end,:)').*cleanNaN(bayesmlpgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlpgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlpgdp(2:end,:)'),2)),2))';
   pooldensforecastpgdp(3,:) = (sum(cleanNaN(bayesdensforecastpgdp(2:end,:)').*cleanNaN(bayesmlpgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlpgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlpgdp(2:end,:)'),2)),2))';
   pooldensforecastpgdp(4,:) = sum(bmamlpgdp.*bmadensforecastpgdp,1);
   
   scorepooldensforecastpgdp(1,:) = (sum(cleanNaN(scoredensforecastpgdp(2:end,:)'),2)./size(cleanNaN(densforecastpgdp(2:k+1,:)'),2))';
   scorepooldensforecastpgdp(2,:) = (sum(cleanNaN(scoredensforecastpgdp(2:end,:)').*cleanNaN(bayesmlpgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlpgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlpgdp(2:end,:)'),2)),2))';
   scorepooldensforecastpgdp(3,:) = (sum(cleanNaN(scorebayesdensforecastpgdp(2:end,:)').*cleanNaN(bayesmlpgdp(2:end,:)')./repmat(sum(cleanNaN(bayesmlpgdp(2:end,:)'),2),1,size(cleanNaN(bayesmlpgdp(2:end,:)'),2)),2))';
   scorepooldensforecastpgdp(4,:) = sum(bmamlpgdp.*scorebmadensforecastpgdp,1);

   
   save(strcat('PITSgdp',num2str(hzn),'.mat'), 'densforecastgdp', 'pooldensforecastgdp', 'factordensforecastgdp', 'vardensforecastgdp',...
                                               'scoredensforecastgdp', 'scorepooldensforecastgdp', 'scorefactordensforecastgdp', 'scorevardensforecastgdp','dategdp', 'datepgdp', 'vnamegdp');
   save(strcat('PITSpgdp',num2str(hzn),'.mat'), 'densforecastpgdp', 'pooldensforecastpgdp', 'factordensforecastpgdp', 'vardensforecastpgdp',...
                                                'scoredensforecastpgdp', 'scorepooldensforecastpgdp', 'scorefactordensforecastpgdp', 'scorevardensforecastpgdp','datepgdp', 'vnamepgdp');
end

disp('Program Complete')