%
%  Matlab file for evaluation of Blanchard-Quah based on the
%  rbc model with stochastic growth
%
%  by: Alain Guay
%      April 2010
  

 zz=1;    %   zz=1  productivity in difference
        %   zz=2  output in difference
        
%  Standard errors of technology shock (sigmav), government spending shocks (sigmag) and
%  preference shocks (sigmap)

sigmav=0.01;
sigmag=0.00;
sigmap=0.01;

nv=2;  % number of variables
nger=250;   %number of observations

% First order correlation fof government spending process (rhog) and preference process (rhop) 

rhop=.95;
rhog=.95;

% number of replications

nsim=50;

% number of impulses for the estimation of the long run effect

nimp=250;

omegabar=2*pi/150;
%omegainf=-omegabar;
omegainf=-omegabar;
omegasup=omegabar;

% Horizon for the impulse responses

horizon=200;
nplot=horizon+1;
nord=0:1:nplot-1;
nord=nord';
hh=12;

% Here, we call the function file to solve the model

[PI,C,Mke]=solvemrbc(sigmav,sigmag,sigmap,rhog,rhop);

[imptrue11,imptrue12,imptrue21,imptrue22,imp11,imp12,MP,CC,PIS]=impptrbc(PI,Mke,sigmav,sigmag,sigmap,zz,nimp);
A0M=[imptrue11(1) imptrue12(1);imptrue21(1) imptrue22(1)];
[vars11,vars12,vars21,vars22]=vardensity(imp11,imp12,imptrue21,imptrue22,nimp,omegasup,omegainf)


implr11=zeros(nimp,nsim); implr12=zeros(nimp,nsim); implr21=zeros(nimp,nsim); implr22=zeros(nimp,nsim);
impfreq11=zeros(nimp,nsim); impfreq12=zeros(nimp,nsim); impfreq21=zeros(nimp,nsim); impfreq22=zeros(nimp,nsim);
impfirst11=zeros(nimp,nsim); impfirst12=zeros(nimp,nsim); impfirst21=zeros(nimp,nsim); impfirst22=zeros(nimp,nsim);
impDO11=zeros(nimp,nsim); impDO12=zeros(nimp,nsim); impDO21=zeros(nimp,nsim); impDO22=zeros(nimp,nsim);



% [sp11,sp12,sp21,sp22,ww]=spectral(imp11,imp12,imptrue21,imptrue22,nimp);
% 
% 
% clf
% subplot(2,2,1),plot(ww,sp11),title('spectral density of shock 1 for variable 1 in level'),
% subplot(2,2,2),plot(ww,sp12),title('spectral density of shock 2 for variable 1 in level'),
% subplot(2,2,3),plot(ww,sp21),title('spectral density of shock 1 for variable 2 in level'),
% subplot(2,2,4),plot(ww,sp22),title('spectral density of shock 2 for variable 2 in level'),
% pause

[sp11,sp12,sp21,sp22,ww]=spectral_window(imp11,imp12,imptrue21,imptrue22,nimp,omegainf,omegasup);

subplot(2,2,1),plot(ww,sp11),title('spectral density of shock 1 for variable 1 in level for the selected window'),
subplot(2,2,2),plot(ww,sp12),title('spectral density of shock 2 for variable 1 in level'),
subplot(2,2,3),plot(ww,sp21),title('spectral density of shock 1 for variable 2 in level'),
subplot(2,2,4),plot(ww,sp22),title('spectral density of shock 2 for variable 2 in level'),
pause


tlr=zeros(nsim,1);tfirst=zeros(nsim,1);tfreq=zeros(nsim,1);jstat=zeros(nsim,1);
corrlr11=zeros(nsim,1);corrfirst11=zeros(nsim,1);corrfreq11=zeros(nsim,1);corrDO11=zeros(nsim,1);
corrlr12=zeros(nsim,1);corrfirst12=zeros(nsim,1);corrfreq12=zeros(nsim,1);corrDO12=zeros(nsim,1);
corrlr21=zeros(nsim,1);corrfirst21=zeros(nsim,1);corrfreq21=zeros(nsim,1);corrDO21=zeros(nsim,1);
corrlr22=zeros(nsim,1);corrfirst22=zeros(nsim,1);corrfreq22=zeros(nsim,1);corrDO22=zeros(nsim,1);

npv5=0;
npv10=0;

nl=4;  % nb of lags in the VAR

rnorm=randn(50000,(nv^2)*nl);
rchi=rnorm.^2;
    
randn('seed',43276);

for i=1:nsim
    i
    [DXN,Nlabor,dlabor,lcy,labor,esimul,PIS,da]=generrbc(PI,Mke,sigmav,sigmag,sigmap,nger,zz);
    dX=[da Nlabor];
  
    [nobs,nv]=size(dX);

    [B,E,sigma,sigbeta,y,y1]=vectar(dX,nl);
    [msim,nnsim]=size(esimul);
    [me,ne]=size(E);
    
    esimm=esimul(msim-me+1:msim,:);
    etech=esimm(:,1);
    etrans=esimm(:,3);
    
    [C,A,D]=ssvar(B,sigma,nv,nl);

    % computation of impulse response

    % Impulse responses from the reduced form
    [IMP]=impulse(C,A,D,nimp);

    % Impulse responses based on long run restrictions
    [A0,A1,IMPA,CIMPA,C1]=bq(IMP,sigma,nv,nimp,A,nl);
    tlr(i)=A0(2,1);
    elr=inv(A0)*E';
    elr=elr';
    %elr=E*A0;
    elr1=elr(:,1);
    elr2=elr(:,2);
    corrlr=corrcoef(elr1,etech);
    corrlr11(i)=corrlr(1,2);
    corrlr=corrcoef(elr1,etrans);
    corrlr12(i)=corrlr(1,2);
    corrlr=corrcoef(elr2,etech);
    corrlr21(i)=corrlr(1,2);
    corrlr=corrcoef(elr2,etrans);
    corrlr22(i)=corrlr(1,2);
    
    A0

    % DiCecio-Owyang estimator
    [VDO,lambda]=spectral_window_reduced(IMP,sigma,nimp,omegainf,omegasup)
    
    [IMPADO,CIMPADO]=impulsefreq(IMP,VDO,nimp,nv);
    eDO=inv(VDO)*E';
    eDO=eDO';
    %eDO=E*VDO;
    eDO1=eDO(:,1);
    eDO2=eDO(:,2);
    corrDO=corrcoef(eDO1,etech);
    corrDO11(i)=corrDO(1,2);
    corrDO=corrcoef(eDO1,etrans);
    corrDO12(i)=corrDO(1,2);
    corrDO=corrcoef(eDO2,etech);
    corrDO21(i)=corrDO(1,2);
    corrDO=corrcoef(eDO2,etrans);
    corrDO22(i)=corrDO(1,2);
    
    % First step estimator
    [theta1]=firststep(IMP,omegasup,omegainf,nimp);
    b=theta1
    [theta]=Afct(A0,b,sigma);
    A0H1=[theta(1) b*theta(2); theta(3) theta(2)]
    tfirst(i)=A0H1(2,1);
    [IMPAfirst,CIMPAfirst]=impulsefreq(IMP,A0H1,nimp,nv);
    efirst=inv(A0H1)*E';
    efirst=efirst';
    %efirst=E*A0H1;
    efirst1=efirst(:,1);
    efirst2=efirst(:,2);
    corrfirst=corrcoef(efirst1,etech);
    corrfirst11(i)=corrfirst(1,2);
    corrfirst=corrcoef(efirst1,etrans);
    corrfirst12(i)=corrfirst(1,2);
    corrfirst=corrcoef(efirst2,etech);
    corrfirst21(i)=corrfirst(1,2);
    corrfirst=corrcoef(efirst2,etrans);
    corrfirst22(i)=corrfirst(1,2);
    
    % Derivatives of impulse responses wrt the VAR parameters
    [derC11,derC12]=deriv(A,nl,nimp,nv);
    %[G,derC11,derC12]=derlutk(A,IMP,nl,nimp,nv);
    
    theta1=b;
    [Cmatrix2,Cmatrix]=coeffals(derC11,derC12,omegasup,omegainf,sigbeta,nimp,theta1); 
    %[Cmatrix2,Cmatrix]=coeffals_id(derC11,derC12,omegasup,omegainf,sigbeta,nimp,theta1);
    [ss1,ss2]=sfct1(IMP,derC11,derC12,omegasup,omegainf,sigbeta,nimp,theta1);
    %[ss1,ss2]=sfct_id(IMP,derC11,derC12,omegasup,omegainf,sigbeta,nimp,theta1)
    [b,fctals,iCmatrix,sf,sds,sdstik]=salstik(Cmatrix2,ss1,ss2,omegasup,(nv^2)*nl);
    betafreq=b
    %end
    [theta]=Afct(A0H1,b,sigma);
    A0H=[theta(1) b*theta(2); theta(3) theta(2)]
    tfreq(i)=A0H(2,1);
    efreq=inv(A0H)*E';
    efreq=efreq';
    %efreq=E*A0H;
    efreq1=efreq(:,1);
    efreq2=efreq(:,2);
    corrfreq=corrcoef(efreq1,etech);
    corrfreq11(i)=corrfreq(1,2);
    corrfreq=corrcoef(efreq1,etrans);
    corrfreq12(i)=corrfreq(1,2);
    corrfreq=corrcoef(efreq2,etech);
    corrfreq21(i)=corrfreq(1,2);
    corrfreq=corrcoef(efreq2,etrans);
    corrfreq22(i)=corrfreq(1,2);
    
    
    %[jstat(i),sds,sdstik,iCmatrixj]=jstatwtik(derC11,derC12,omegasup,omegainf,sigbeta,nimp,betafreq,nobs,sf,iCmatrix);
    %jstat(i) = nobs*sf'*iCmatrix*sf;
    vv=sds./sdstik;
    %vv(1)=0;
    rv=rchi*vv;
    srv=sort(rv);
    cvtik5=srv(.95*50000);
    cvtik10=srv(.90*50000);
    jstat(i)=nobs*fctals;
    jstati=jstat(i)
    if jstat(i) >= cvtik5
       npv5=npv5+1;
    end
    if jstat(i) >= cvtik10
       npv10=npv10+1;
    end   
    [IMPAfreq,CIMPAfreq]=impulsefreq(IMP,A0H,nimp,nv);
    %gbar=sf;
    %[Cmatrix2s]=coeffals(derC11,derC12,omegasup,omegainf,sigbeta,nimp,b);
    
      
    for j=1:nv
      numj=int2str(j);
      for k=1:nv
        numk=int2str(k);
        for ll=0:nimp
            eval(['impfreq',numj,numk,'(ll+1,i)=IMPAfreq(j,k,ll+1);'])
            eval(['implr',numj,numk,'(ll+1,i)=IMPA(j,k,ll+1);'])
            eval(['impfirst',numj,numk,'(ll+1,i)=IMPAfirst(j,k,ll+1);'])
            eval(['impDO',numj,numk,'(ll+1,i)=IMPADO(j,k,ll+1);'])
        end
        if j==1 & k==1 
           impfreq11(:,i)=cumsum(impfreq11(:,i));
           implr11(:,i)=cumsum(implr11(:,i));
           impfirst11(:,i)=cumsum(impfirst11(:,i));
           impDO11(:,i)=cumsum(impDO11(:,i));
           
        elseif j==1 & k==2   
           impfreq12(:,i)=cumsum(impfreq12(:,i));
           implr12(:,i)=cumsum(implr12(:,i));
           impfirst12(:,i)=cumsum(impfirst12(:,i));
           impDO12(:,i)=cumsum(impDO12(:,i));
        end
       %eval(['implfreq=impfreq',numj,numk,';'])
       %eval(['implr=implr',numj,numk,';'])
      % plot(nord,implfreq(1:nplot) implt(1:nplot)]),title(['Freq impulse response of shock ', num2str(k),'  on variable ', num2str(j),' in level'])
       %pause
     end
   end

%     plot(nord,[imptrue11(1:nplot) implr11(1:nplot,i) impfreq11(1:nplot,i)] ),title(['impulse response of shock 1 on variable 1 in level'])
%     pause        
% 
%     plot(nord,[imptrue12(1:nplot) implr12(1:nplot,i) impfreq12(1:nplot,i)] ),title(['impulse response of shock 2 on variable 1 in level'])
%     pause        
% 
%     plot(nord,[imptrue21(1:nplot) implr21(1:nplot,i) impfreq21(1:nplot,i)] ),title(['impulse response of shock 1 on variable 2 in level'])
%     pause        
% 
%     plot(nord,[imptrue22(1:nplot) implr22(1:nplot,i) impfreq22(1:nplot,i)] ),title(['impulse response of shock 2 on variable 2 in level'])
%     pause        
    
end

%     for j=1:nv
%       numj=int2str(j);
%       for k=1:nv
%         numk=int2str(k);
%         for ll=0:nimp
%             eval(['impfreq',numj,numk,'(ll+1)=IMPA(j,k,ll+1);'])
%         end
%         if j==1 & k==1 
%            impfreq11=cumsum(impfreq11);
%            
%         elseif j==1 & k==2   
%            impfreq12=cumsum(impfreq12);
%         end
%        eval(['implfreq=impfreq',numj,numk,';'])
%        plot(nord,implfreq(1:nplot)),title(['Freq impulse response of shock ', num2str(k),'  on variable ', num2str(j),' in level'])
%        pause
%      end
%    end

%     plot(nord,[imptrue11(1:nplot) implr11(1:nplot)' impfreq11(1:nplot)'] ),title(['impulse response of shock 1 on variable 1 in level'])
%     pause        
% 
%     plot(nord,[imptrue12(1:nplot) implr12(1:nplot)' impfreq12(1:nplot)'] ),title(['impulse response of shock 2 on variable 1 in level'])
%     pause        
% 
%     plot(nord,[imptrue21(1:nplot) implr21(1:nplot)' impfreq21(1:nplot)'] ),title(['impulse response of shock 1 on variable 2 in level'])
%     pause        
% 
%     plot(nord,[imptrue22(1:nplot) implr22(1:nplot)' impfreq22(1:nplot)'] ),title(['impulse response of shock 2 on variable 2 in level'])
%     pause

npv5=npv5/nsim
npv10=npv10/nsim
mfreq=mean(tfreq);
mfirst=mean(tfirst);
mlr=mean(tlr);
d=imptrue21(1);
msefreq=(mfreq-d)^2+var(tfreq)
msefirst=(mfirst-d)^2+var(tfirst)
mselr=(mlr-d)^2+var(tlr)

mcorrlr11=mean(corrlr11);
mcorrfirst11=mean(corrfirst11);
mcorrfreq11=mean(corrfreq11);
mcorrDO11=mean(corrDO11);
mcorrlr12=mean(corrlr12);
mcorrfirst12=mean(corrfirst12);
mcorrfreq12=mean(corrfreq12);
mcorrDO12=mean(corrDO12);
mcorrlr21=mean(corrlr21);
mcorrfirst21=mean(corrfirst21);
mcorrfreq21=mean(corrfreq21);
mcorrDO21=mean(corrDO21);
mcorrlr22=mean(corrlr22);
mcorrfirst22=mean(corrfirst22);
mcorrfreq22=mean(corrfreq22);
mcorrDO22=mean(corrDO22);


mimplr11=mean(implr11')';
simplr11=sort(implr11');
simplr115=simplr11(ceil(.05*nsim),:);
simplr1195=simplr11(ceil(.95*nsim),:);

erlr11=implr11(1:hh,:)-imptrue11(1:hh)*ones(1,nsim);
serlr11=erlr11.^2;
merlr11=mean((erlr11'));
merlr11=abs(merlr11);
smerlr11=[merlr11(1) sum(merlr11(1:5)) sum(merlr11(1:9)) sum(merlr11(1:hh))];
rmselr11=sqrt(merlr11.^2+var(implr11(1:hh,:)'));
%rmselr11=sqrt(mean(serlr11));
srmselr11=[rmselr11(1) sum(rmselr11(1:5)) sum(rmselr11(1:9)) sum(rmselr11(1:hh))];

clf
plot(nord,[imptrue11(1:nplot) mimplr11(1:nplot) simplr115(1:nplot)' simplr1195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 1 in level: LR'])
    pause        

mimpfirst11=mean(impfirst11')';
simpfirst11=sort(impfirst11');
simpfirst115=simpfirst11(ceil(.05*nsim),:);
simpfirst1195=simpfirst11(ceil(.95*nsim),:);

erfirst11=impfirst11(1:hh,:)-imptrue11(1:hh)*ones(1,nsim);
serfirst11=erfirst11.^2;
merfirst11=mean((erfirst11'));
merfirst11=abs(merfirst11);
smerfirst11=[merfirst11(1) sum(merfirst11(1:5)) sum(merfirst11(1:9)) sum(merfirst11(1:hh))];
rmsefirst11=sqrt(merfirst11.^2+var(impfirst11(1:hh,:)'));
%rmsefirst11=sqrt(mean(serfirst11));
srmsefirst11=[rmsefirst11(1) sum(rmsefirst11(1:5)) sum(rmsefirst11(1:9)) sum(rmsefirst11(1:hh))];


plot(nord,[imptrue11(1:nplot) mimpfirst11(1:nplot) simpfirst115(1:nplot)' simpfirst1195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 1 in level: First'])
    pause       
    
mimpfreq11=mean(impfreq11')';
simpfreq11=sort(impfreq11');
simpfreq115=simpfreq11(ceil(.05*nsim),:);
simpfreq1195=simpfreq11(ceil(.95*nsim),:);

erfreq11=impfreq11(1:hh,:)-imptrue11(1:hh)*ones(1,nsim);
serfreq11=erfreq11.^2;
merfreq11=mean((erfreq11'));
merfreq11=abs(merfreq11);
smerfreq11=[merfreq11(1) sum(merfreq11(1:5)) sum(merfreq11(1:9)) sum(merfreq11(1:hh))];
rmsefreq11=sqrt(merfreq11.^2+var(impfreq11(1:hh,:)'));
%rmsefreq11=sqrt(mean(serfreq11));
srmsefreq11=[rmsefreq11(1) sum(rmsefreq11(1:5)) sum(rmsefreq11(1:9)) sum(rmsefreq11(1:hh))];


plot(nord,[imptrue11(1:nplot) mimpfreq11(1:nplot) simpfreq115(1:nplot)' simpfreq1195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 1 in level: Freq'])
    pause   
    
mimpDO11=mean(impDO11')';
simpDO11=sort(impDO11');
simpDO115=simpDO11(ceil(.05*nsim),:);
simpDO1195=simpDO11(ceil(.95*nsim),:);

erDO11=impDO11(1:hh,:)-imptrue11(1:hh)*ones(1,nsim);
serDO11=erDO11.^2;
merDO11=mean((erDO11'));
merDO11=abs(merDO11);
smerDO11=[merDO11(1) sum(merDO11(1:5)) sum(merDO11(1:9)) sum(merDO11(1:hh))];
rmseDO11=sqrt(merDO11.^2+var(impDO11(1:hh,:)'));
%rmseDO11=sqrt(mean(serDO11));
srmseDO11=[rmseDO11(1) sum(rmseDO11(1:5)) sum(rmseDO11(1:9)) sum(rmseDO11(1:hh))];


plot(nord,[imptrue11(1:nplot) mimpDO11(1:nplot) simpDO115(1:nplot)' simpDO1195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 1 in level: DO'])
    pause       

mimplr12=mean(implr12')';
simplr12=sort(implr12');
simplr125=simplr12(ceil(.05*nsim),:);
simplr1295=simplr12(ceil(.95*nsim),:);

erlr12=implr12(1:hh,:)-imptrue12(1:hh)*ones(1,nsim);
serlr12=erlr12.^2;
merlr12=mean((erlr12'));
merlr12=abs(merlr12);
smerlr12=[merlr12(1) sum(merlr12(1:5)) sum(merlr12(1:9)) sum(merlr12(1:hh))];
rmselr12=sqrt(merlr12.^2+var(implr12(1:hh,:)'));
%rmselr12=sqrt(mean(serlr12));
srmselr12=[rmselr12(1) sum(rmselr12(1:5)) sum(rmselr12(1:9)) sum(rmselr12(1:hh))];

plot(nord,[imptrue12(1:nplot) mimplr12(1:nplot) simplr125(1:nplot)' simplr1295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 1 in level: LR'])
    pause        

mimpfirst12=mean(impfirst12')';
simpfirst12=sort(impfirst12');
simpfirst125=simpfirst12(ceil(.05*nsim),:);
simpfirst1295=simpfirst12(ceil(.95*nsim),:);

erfirst12=impfirst12(1:hh,:)-imptrue12(1:hh)*ones(1,nsim);
serfirst12=erfirst12.^2;
merfirst12=mean((erfirst12'));
merfirst12=abs(merfirst12);
smerfirst12=[merfirst12(1) sum(merfirst12(1:5)) sum(merfirst12(1:9)) sum(merfirst12(1:hh))];
rmsefirst12=sqrt(merfirst12.^2+var(impfirst12(1:hh,:)'));
%rmsefirst12=sqrt(mean(serfirst12));
srmsefirst12=[rmsefirst12(1) sum(rmsefirst12(1:5)) sum(rmsefirst12(1:9)) sum(rmsefirst12(1:hh))];

plot(nord,[imptrue12(1:nplot) mimpfirst12(1:nplot) simpfirst125(1:nplot)' simpfirst1295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 1 in level: First'])
    pause       
    
mimpfreq12=mean(impfreq12')';
simpfreq12=sort(impfreq12');
simpfreq125=simpfreq12(ceil(.05*nsim),:);
simpfreq1295=simpfreq12(ceil(.95*nsim),:);

erfreq12=impfreq12(1:hh,:)-imptrue12(1:hh)*ones(1,nsim);
serfreq12=erfreq12.^2;
merfreq12=mean((erfreq12'));
merfreq12=abs(merfreq12);
smerfreq12=[merfreq12(1) sum(merfreq12(1:5)) sum(merfreq12(1:9)) sum(merfreq12(1:hh))];
rmsefreq12=sqrt(merfreq12.^2+var(impfreq12(1:hh,:)'));
%rmsefreq11=sqrt(mean(serfreq11));
srmsefreq12=[rmsefreq12(1) sum(rmsefreq12(1:5)) sum(rmsefreq12(1:9)) sum(rmsefreq12(1:hh))];


plot(nord,[imptrue12(1:nplot) mimpfreq12(1:nplot) simpfreq125(1:nplot)' simpfreq1295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 1 in level: Freq'])
    pause   
    
    
mimpDO12=mean(impDO12')';
simpDO12=sort(impDO12');
simpDO125=simpDO12(ceil(.05*nsim),:);
simpDO1295=simpDO12(ceil(.95*nsim),:);

erDO12=impDO12(1:hh,:)-imptrue12(1:hh)*ones(1,nsim);
serDO12=erDO12.^2;
merDO12=mean((erDO12'));
merDO12=abs(merDO12);
smerDO12=[merDO12(1) sum(merDO12(1:5)) sum(merDO12(1:9)) sum(merDO12(1:hh))];
rmseDO12=sqrt(merDO12.^2+var(impDO12(1:hh,:)'));
%rmseDO12=sqrt(mean(serDO12));
srmseDO12=[rmseDO12(1) sum(rmseDO12(1:5)) sum(rmseDO12(1:9)) sum(rmseDO12(1:hh))];

plot(nord,[imptrue12(1:nplot) mimpDO12(1:nplot) simpDO125(1:nplot)' simpDO1295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 1 in level DO'])
    pause   

mimplr21=mean(implr21')';
simplr21=sort(implr21');
simplr215=simplr21(ceil(.05*nsim),:);
simplr2195=simplr21(ceil(.95*nsim),:);

erlr21=implr21(1:hh,:)-imptrue21(1:hh)*ones(1,nsim);
serlr21=erlr21.^2;
merlr21=mean((erlr21)');
merlr21=abs(merlr21);
smerlr21=[merlr21(1) sum(merlr21(1:5)) sum(merlr21(1:9)) sum(merlr21(1:hh))];
rmselr21=sqrt(merlr21.^2+var(implr21(1:hh,:)'));
%rmselr21=sqrt(mean(serlr21));
srmselr21=[rmselr21(1) sum(rmselr21(1:5)) sum(rmselr21(1:9)) sum(rmselr21(1:hh))];

plot(nord,[imptrue21(1:nplot) mimplr21(1:nplot) simplr215(1:nplot)' simplr2195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 2 in level: LR'])
    pause        

mimpfirst21=mean(impfirst21')';
simpfirst21=sort(impfirst21');
simpfirst215=simpfirst21(ceil(.05*nsim),:);
simpfirst2195=simpfirst21(ceil(.95*nsim),:);

erfirst21=impfirst21(1:hh,:)-imptrue21(1:hh)*ones(1,nsim);
serfirst21=erfirst21.^2;
merfirst21=mean((erfirst21'));
merfirst21=abs(merfirst21);
smerfirst21=[merfirst21(1) sum(merfirst21(1:5)) sum(merfirst21(1:9)) sum(merfirst21(1:hh))];
rmsefirst21=sqrt(merfirst21.^2+var(impfirst21(1:hh,:)'));
%sqrt(mean(serfirst21));
srmsefirst21=[rmsefirst21(1) sum(rmsefirst21(1:5)) sum(rmsefirst21(1:9)) sum(rmsefirst21(1:hh))];

plot(nord,[imptrue21(1:nplot) mimpfirst21(1:nplot) simpfirst215(1:nplot)' simpfirst2195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 2 in level: First'])
    pause      
    
mimpfreq21=mean(impfreq21')';
simpfreq21=sort(impfreq21');
simpfreq215=simpfreq21(ceil(.05*nsim),:);
simpfreq2195=simpfreq21(ceil(.95*nsim),:);

erfreq21=impfreq21(1:hh,:)-imptrue21(1:hh)*ones(1,nsim);
serfreq21=erfreq21.^2;
merfreq21=mean((erfreq21'));
merfreq21=abs(merfreq21);
smerfreq21=[merfreq21(1) sum(merfreq21(1:5)) sum(merfreq21(1:9)) sum(merfreq21(1:hh))];
rmsefreq21=sqrt(merfreq21.^2+var(impfreq21(1:hh,:)'));
%sqrt(mean(serfreq21));
srmsefreq21=[rmsefreq21(1) sum(rmsefreq21(1:5)) sum(rmsefreq21(1:9)) sum(rmsefreq21(1:hh))];

plot(nord,[imptrue21(1:nplot) mimpfreq21(1:nplot) simpfreq215(1:nplot)' simpfreq2195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 2 in level: Freq'])
    pause   

mimpDO21=mean(impDO21')';
simpDO21=sort(impDO21');
simpDO215=simpDO21(ceil(.05*nsim),:);
simpDO2195=simpDO21(ceil(.95*nsim),:);

erDO21=impDO21(1:hh,:)-imptrue21(1:hh)*ones(1,nsim);
serDO21=erDO21.^2;
merDO21=mean((erDO21'));
merDO21=abs(merDO21);
smerDO21=[merDO21(1) sum(merDO21(1:5)) sum(merDO21(1:9)) sum(merDO21(1:hh))];
rmseDO21=sqrt(merDO21.^2+var(impDO21(1:hh,:)'));
%sqrt(mean(serDO21));
srmseDO21=[rmseDO21(1) sum(rmseDO21(1:5)) sum(rmseDO21(1:9)) sum(rmseDO21(1:hh))];

plot(nord,[imptrue21(1:nplot) mimpDO21(1:nplot) simpDO215(1:nplot)' simpDO2195(1:nplot)' ] ),title(['impulse response of shock 1 on variable 2 in level: DO'])
    pause   

    
mimplr22=mean(implr22')';
simplr22=sort(implr22');
simplr225=simplr22(ceil(.05*nsim),:);
simplr2295=simplr22(ceil(.95*nsim),:);

erlr22=implr22(1:hh,:)-imptrue22(1:hh)*ones(1,nsim);
serlr22=erlr22.^2;
merlr22=mean((erlr22'));
merlr22=abs(merlr22);
smerlr22=[merlr22(1) sum(merlr22(1:5)) sum(merlr22(1:9)) sum(merlr22(1:hh))];
rmselr22=sqrt(merlr22.^2+var(implr22(1:hh,:)'));
%rmselr22=sqrt(mean(serlr22));
srmselr22=[rmselr22(1) sum(rmselr22(1:5)) sum(rmselr22(1:9)) sum(rmselr22(1:hh))];


plot(nord,[imptrue22(1:nplot) mimplr22(1:nplot) simplr225(1:nplot)' simplr2295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 2 in level: LR'])
    pause  
    
mimpfirst22=mean(impfirst22')';
simpfirst22=sort(impfirst22');
simpfirst225=simpfirst22(ceil(.05*nsim),:);
simpfirst2295=simpfirst22(ceil(.95*nsim),:);

erfirst22=impfirst22(1:hh,:)-imptrue22(1:hh)*ones(1,nsim);
serfirst22=erfirst22.^2;
merfirst22=mean((erfirst22'));
merfirst22=abs(merfirst22);
smerfirst22=[merfirst22(1) sum(merfirst22(1:5)) sum(merfirst22(1:9)) sum(merfirst22(1:hh))];
rmsefirst22=sqrt(merfirst22.^2+var(impfirst22(1:hh,:)'));
%rmsefirst22=sqrt(mean(serfirst22));
srmsefirst22=[rmsefirst22(1) sum(rmsefirst22(1:5)) sum(rmsefirst22(1:9)) sum(rmsefirst22(1:hh))];

plot(nord,[imptrue22(1:nplot) mimpfirst22(1:nplot) simpfirst225(1:nplot)' simpfirst2295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 2 in level: First'])
    pause   
    

mimpfreq22=mean(impfreq22')';
simpfreq22=sort(impfreq22');
simpfreq225=simpfreq22(ceil(.05*nsim),:);
simpfreq2295=simpfreq22(ceil(.95*nsim),:);

erfreq22=impfreq22(1:hh,:)-imptrue22(1:hh)*ones(1,nsim);
serfreq22=erfreq22.^2;
merfreq22=mean((erfreq22'));
merfreq22=abs(merfreq22);
smerfreq22=[merfreq22(1) sum(merfreq22(1:5)) sum(merfreq22(1:9)) sum(merfreq22(1:hh))];
rmsefreq22=sqrt(merfreq22.^2+var(impfreq22(1:hh,:)'));
%rmsefreq22=sqrt(mean(serfreq22));
srmsefreq22=[rmsefreq22(1) sum(rmsefreq22(1:5)) sum(rmsefreq22(1:9)) sum(rmsefreq22(1:hh))];

plot(nord,[imptrue22(1:nplot) mimpfreq22(1:nplot) simpfreq225(1:nplot)' simpfreq2295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 2 in level: Freq'])
    pause   
    
mimpDO22=mean(impDO22')';
simpDO22=sort(impDO22');
simpDO225=simpDO22(ceil(.05*nsim),:);
simpDO2295=simpDO22(ceil(.95*nsim),:);

erDO22=impDO22(1:hh,:)-imptrue22(1:hh)*ones(1,nsim);
serDO22=erDO22.^2;
merDO22=mean((erDO22'));
merDO22=abs(merDO22);
smerDO22=[merDO22(1) sum(merDO22(1:5)) sum(merDO22(1:9)) sum(merDO22(1:hh))];
rmseDO22=sqrt(merDO22.^2+var(impDO22(1:hh,:)'));
%rmseDO22=sqrt(mean(serDO22));
srmseDO22=[rmseDO22(1) sum(rmseDO22(1:5)) sum(rmseDO22(1:9)) sum(rmseDO22(1:hh))];

plot(nord,[imptrue22(1:nplot) mimpDO22(1:nplot) simpDO225(1:nplot)' simpDO2295(1:nplot)' ] ),title(['impulse response of shock 2 on variable 2 in level: DO'])
    pause     
    
disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  Correlation of the first shock with the true permanent shock  ');
disp('________________________________________________________________');
disp('   LR        FIRST     FREQUENCY   DO ');
disp('________________________________________________________________');
disp('');
disp( [mcorrlr11 mcorrfirst11 mcorrfreq11  mcorrDO11] );
disp('________________________________________________________________');   

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  Correlation of the first shock with the true transitory shock  ');
disp('________________________________________________________________');
disp('   LR        FIRST     FREQUENCY  DO');
disp('________________________________________________________________');
disp('');
disp( [mcorrlr12 mcorrfirst12 mcorrfreq12 mcorrDO12] );
disp('________________________________________________________________');    

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  Correlation of the second shock with the true permanent shock  ');
disp('________________________________________________________________');
disp('   LR        FIRST     FREQUENCY  DO');
disp('________________________________________________________________');
disp('');
disp( [mcorrlr21 mcorrfirst21 mcorrfreq21 mcorrDO21] );
disp('________________________________________________________________');    

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  Correlation of the second shock with the true transitory shock  ');
disp('________________________________________________________________');
disp('   LR        FIRST     FREQUENCY  DO');
disp('________________________________________________________________');
disp('');
disp( [mcorrlr22 mcorrfirst22 mcorrfreq22 mcorrDO22] );
disp('________________________________________________________________');    

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: implr11  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12  ');
disp('________________________________________________________________');
disp('');
disp( smerlr11);
disp( srmselr11);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfirst11  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfirst11);
disp( srmsefirst11);
disp('________________________________________________________________');


disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfreq11  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfreq11);
disp( srmsefreq11);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impDO11  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerDO11);
disp( srmseDO11);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: implr12  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12  ');
disp('________________________________________________________________');
disp('');
disp( smerlr12);
disp( srmselr12);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean square errors: impfirst12  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfirst12);
disp( srmsefirst12);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfreq12  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfreq12);
disp( srmsefreq12);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impDO12  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerDO12);
disp( srmseDO12);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: implr21  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12  ');
disp('________________________________________________________________');
disp('');
disp( smerlr21);
disp( srmselr21);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfirst21  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfirst21);
disp( srmsefirst21);
disp('________________________________________________________________');


disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfreq21  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfreq21);
disp( srmsefreq21);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impDO21  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerDO21);
disp( srmseDO21);
disp('________________________________________________________________');



disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: implr22  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12  ');
disp('________________________________________________________________');
disp('');
disp( smerlr22);
disp( srmselr22);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfirst22  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfirst22);
disp( srmsefirst22);
disp('________________________________________________________________');


disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impfreq22  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerfreq22);
disp( srmsefreq22);
disp('________________________________________________________________');

disp('----------------------------------------------------------------');
disp('');
disp('________________________________________________________________');
disp('  mean absolute bias and root mean squre errors: impDO22  ');
disp('________________________________________________________________');
disp('   0          0 to 4     0 to 8     0 to 12 ');
disp('________________________________________________________________');
disp('');
disp( smerDO22);
disp( srmseDO22);
disp('________________________________________________________________');