%ROLLING NESTED
%This runs MC simulations for the null hypotheses in the paper -- select
%the null hypotheses by commenting out the ones you do not want in lines 60-62
clear; randn('seed',13); 
MC=5000; 5000; 5000;
table=[]; table1=[];  h=1; table2=[];

DGP=2; %=1 if serial correlation and =2 if iid
wantaltern=2; %wantaltern=1 if want size, =2 if TVP, =3 if pred content, =4 if overfit

if wantaltern==1; Rvec=[20,50,100,200];  m=100; Predv=[m+50, m+100, m+200];  alternativeV=0; end;   
if wantaltern>1;  Rvec=100; Predv=300; m=100;  end;
if wantaltern==2; m=60; if DGP==2; alternativeV=[0:0.1:1]; elseif DGP==1; alternativeV=[0:0.1:1]; end; end;
if wantaltern==3; alternativeV=[0:0.1:1]; end; 
if wantaltern==4; alternativeV=[0,1,2,5:5:40]; end; 

if DGP==1; 
    str=['delete main_size_rollDGP1alternative',num2str(wantaltern),'.out; diary main_size_rollDGP1alternative',num2str(wantaltern),'.out;']; 
elseif DGP==2; 
    str=['delete main_size_rollDGP2alternative',num2str(wantaltern),'.out; diary main_size_rollDGP2alternative',num2str(wantaltern),'.out;']; 
end;
eval(str);

for nullhp=[1,2,3];
    for alternative=alternativeV; 

        for R=Rvec;

            for Pred=Predv;  
                rejv=[]; corrv=[]; betav=[];
                T=R+Pred+1; 

                % Model:
                % y(t)=a2*x(t-1)+e(t)
                % x(t)=b1+b2*x(t-1)+b3*x(t-2)+b4*x(t-3)+b5*x(t-4)+b6*x(t-5)+v(t)
                % e(t)=rho*e(t-1)+eps(t)
                % a2= 0.1461
                % rho=-0.0073
                % var(eps)=2.6363   
                % var(v)=1.5151 not used

                if DGP==1; 
                    b1=0.1409;   
                    b2=-0.1158;
                    b3=0.1059;
                    b4=0.0957;
                    b5=0.0089;
                    b6=0.1412; %modulus(roots([1, -b2, -b3, -b4, -b5, -b6]))
                    b=[b2;b3;b4;b5;b6]; 
                    rho=-0.0073;
                    sigmaeps=2.6363; 
                    sigmae=sigmaeps/(1-rho^2); 
                    qn=2;
                elseif DGP==2; 
                    b1=0; 0.1409;   
                    b2=0; -0.1158;
                    b3=0; 0.1059;
                    b4=0; 0.0957;
                    b5=0; 0.0089;
                    b6=0; 0.1412; %modulus(roots([1, -b2, -b3, -b4, -b5, -b6]))
                    b=[b2;b3;b4;b5;b6]; 
                    rho=0; -0.0073;
                    sigmaeps=1; 2.6363; 
                    sigmae=sigmaeps/(1-rho^2); 
                    qn=0;
                end;

                trash=100; 
                alpha=ABDcalcNEW(R,Pred,DGP,MC); 
                if nullhp==1; 
                a2=alpha(1,1); %H0 for A
                elseif nullhp==2; 
                a2=alpha(1,2); %H0 for B
                elseif nullhp==3; 
                a2=alpha(1,3); %H0 for U
                end;
                
                if wantaltern==2; 
                    test = [1:1:T+1]';
                    a2 = a2 + alternative*cos(3.14*1.5*test./T).*(1-test./T); 
                end;
                if wantaltern==3; a2=a2+alternative; end;

            for rep=1:MC; 
                if DGP==1; x1=recserar(randn(T+1+trash,1)+b1,zeros(1,1),b);  
                elseif DGP==2; x1=ones(T+1+trash,1); end;

                x1=x1(trash+1:end,:); 
                %CONSTANT CASE a2=sqrt(sqrt(3)*(sigmae/(sum(x1(1:R,:).^2)))); %H0 for BETA and for B       (3/(R^2))^(1/4); sqrt(sigma2/sum(x1(1:R).^2));  
                %CONSTANT CASE a2=(sqrt(sigmae)/R)*sqrt(1+sqrt(1+R^2-3*R)); %H0 for U      
                %CONSTANT CASE a2=sqrt(1/R);% sqrt(sigmae/(sum(x1(1:R,:).^2))); %H0 for A      

                e=recserar(randn(T+1+trash,1)*sqrt(sigmae),0,rho); e=e(trash+1:end,:);   
                y=x1.*a2+e;   
                if wantaltern==4; 
                    if DGP==1;  
                    x1T=recserar(randn(T+1+trash,1)+b1,zeros(1,1),b);  
                    elseif DGP==2; 
                    x1T=ones(T+1+trash,1); 
                    end;
                    x1T=x1T(trash+1:end,:); 
                    if alternative==0; x1=x1T;
                    else x1=[x1T,randn(T+1,alternative)]; end; 
                    y=x1T*a2+e; 
                end; 

                %OOS forecasts
                f_in_roll1=[];   f_oos_roll1=[];  b_roll1_v=[];  yf_roll1=[];  true=[];
                f_in_roll2=[];   f_oos_roll2=[];  b_roll2_v=[];  yf_roll2=[];   
                dDeltaL_in=[]; dDeltaL_oos=[]; e_in=[]; f_oos_roll1_true=[];
                for t=R:T; 
                   b_roll1   =olsbeta(  y(t-R+1:t,1),x1(t-R+1:t,:)  );  
                   b_roll2   =0; 
                   true     =[true;    y(t+h,1)];  
                   in1=(( olsres( y(t-R+1:t,1),x1(t-R+1:t,:)  ) ).^2);
                   in2=((  (  y(t-R+1:t,1)  ) ).^2);
                   f_in_roll1    = [f_in_roll1;  in1(end) ]; 
                   f_in_roll2    = [f_in_roll2;  in2(end) ]; 
                   f_oos_roll1   = [f_oos_roll1 ;(y(t+h,1)-x1(t+h,:)  *b_roll1')       .^2]; 
                   f_oos_roll1_true   = [f_oos_roll1_true ;(y(t+h,1)-x1(t+h,:)  *a2)       .^2]; 
                   f_oos_roll2   = [f_oos_roll2 ;(y(t+h,1))       .^2];
                   dDeltaL_in_here = -2*olsres( y(t-R+1:t,1),x1(t-R+1:t,:)  );
                   dDeltaL_in   = [dDeltaL_in; dDeltaL_in_here(end,1)]; 
                   dDeltaL_oos_here = -2*(y(t+h,1)-x1(t+h,:)  *b_roll1');
                   dDeltaL_oos   = [dDeltaL_oos; dDeltaL_oos_here(end,1)]; 
                   e_in_here = olsres( y(t-R+1:t,1),x1(t-R+1:t,:)  );
                   e_in          = [e_in; e_in_here(end,1)]; 
                end; 

                %Our procedure
                DeltaL_oos=f_oos_roll1-f_oos_roll2;     
                DeltaL_in=f_in_roll1-f_in_roll2; 
                results = abc(DeltaL_oos,DeltaL_in,m,qn,dDeltaL_oos,dDeltaL_in,R,e_in); 
                A=results.A; B=results.B; C=results.C; rej=results.rej; 
                rejv=[rejv; rej ]; %betav=[betav; results.beta, var(f_oos_roll1-f_oos_roll2), var(f_oos_roll1_true-f_oos_roll2)];   
            end;  

            resulttable1=[mean(A),mean(rejv(:,1)),mean(B),mean(rejv(:,2)),mean(C),mean(rejv(:,3))]; 
            %resulttable1=[mean(rejv(:,1)),mean(rejv(:,2)),mean(rejv(:,3))]; 
            table1=[table1; alternative,R,Pred,m,resulttable1],
          end; 
        end;
    end;
disp(table1),
            if nullhp==1;     table2=[table2, table1(:,[2:6])]; nrows=size(table1,1); 
            elseif nullhp==2; table2=[table2, table1( nrows+1:end   , [7:8]  )];
            elseif nullhp==3; table2=[table2, table1( 2*nrows+1:end , [9:10] )];
            end;
end;
disp(table2);
diary off; 