%=========
data_all_post1947
dataYY=datahelp(2:end,1);
dataXX=datahelp(1:(end-1),2:16);
size(dataYY)
size(dataXX)
plithosfactors=size(dataXX,2);
plithosmodelon=plithosfactors+1;
holdoutperiod=40; % from 1955:1 to 1964:4 (40 periods)
outofsampleperiod=184; % from 1965:1 to 2010:4 (184 periods)
allperiod=holdoutperiod+outofsampleperiod;


%% MATRICES WITH THE INDIVIDUAL FORECASTS
%==========================================================================
yhat_HO=zeros(allperiod,plithosmodelon); 


%% COMPUTE INDIVIDUAL LINEAR AND QUANTILE REGRESSION FORECATSTS
%==========================================================================
%============================================
%START ITERATIVE FOR THE OUT-OF-SAMPLE PERIOD
%============================================
for iter=1:allperiod
    disp(['Iteration : ' num2str(iter) ' out of ' num2str(allperiod)])
    %============================================================
    %=========
    %Make data
    %=========
    data=[]; datafactors=[];
    data=dataYY(1:(end-allperiod+iter-1));
    datafactors=dataXX(1:(end-allperiod+iter-1),:); 
    datafactors_pred=dataXX(end-allperiod+iter,:);
    [T,plithosmet]=size(data);
    %============================================================
    
    %==== MEAN REGRESSION  ============== 
    %===================================
    %Start iterative work for each model
    %===================================    
    for i=1:plithosmodelon
        model=i-1;     %i=1 refers to a model with constant only
        model(model==0)=[];
        factors_used=[];  factors_used_const=[];
        factors_used=datafactors(:,model);
        numfactors=length(model);
        X=[ones(T,1) factors_used];
        
        X_pred=[];
        X_pred=[1 datafactors_pred(model)];
        
        %calculate conditional mean model
        b=(inv(X'*X))*X'*data;  yhat=X_pred*b;  yhat_HO(iter,i)=yhat;
        
    end
    %=================================
    %End iterative work for each model
    %=================================
end
%==========================================
%END ITERATIVE FOR THE OUT-OF-SAMPLE PERIOD
%==========================================




    if (i/5000)==fix(i/5000), i,  end
    model=all_models(i,:);
    model(model==0)=[];
    
    factors_used=[];  
    factors_used=datafactors(:,model);
    numfactors=length(model);
    
    speck=garchset('P',1,'Q',1,'Regress',[1:numfactors],'display','off');
    Coeff=[]; Errors=[]; LLF=[]; Innovations=[]; Sigmas=[]; Summary=[];
    [Coeff,Errors,LLF,Innovations,Sigmas,Summary]=garchfit(speck,data,factors_used); 
    
    %Calculate AIC, BIC
    loglikelihood(i)=LLF;
    numberparam=1+numfactors+3;   %constant, factor betas, a,b,g of variance equation
    AIC1(i)=-2*LLF+2*(numberparam);
    BIC1(i)=-2*LLF+(log(T)*numberparam);
    
    %Calculate Log-Marginal Likelihood 
    param=[]; covmat=[]; loglik=[]; logprior=[]; 
    %Keep the maximum log-likelihood
    loglik=LLF;

    %Keep the parameters
    param(1)=garchget(Coeff,'C');
    param((1+1):(1+numfactors))=garchget(Coeff,'Regress');
    param(1+numfactors+1)=garchget(Coeff,'K');
    param(1+numfactors+2)=garchget(Coeff,'GARCH');
    param(1+numfactors+3)=garchget(Coeff,'ARCH');
    param=param';
    
    %Keep the covariance matrix BASED ON MLE  (1)
    covmat=Summary.covMatrix;
    stderrors=sqrt(diag(covmat));

    %calculate the log prior
    %Use multivariate Normal prior for the constant m and the betas b1,b2,...,bk
    parbetas=param(1:(1+numfactors));
    mpb=[zeros((1+numfactors),1)];  % mean_prior_betas
    cpb=T*covmat((1:(1+numfactors)),(1:(1+numfactors))); %cov_prior_betas
    logpriorb=[];
    logpriorb=-0.5*(1+numfactors)*log(2*pi)-0.5*log(det(cpb))-0.5*(parbetas-mpb)'*inv(cpb)*(parbetas-mpb);