@ ------------------------------------------- CoRRAM_3.g ---------------------- 

   Alfred Maußner
   28 December 2010

   Third piece of code to solve a DSGE model via linear or quadratic
   approximation of the model at the stationary equilibrium of the
   underlying deterministic model.

----------------------------------------------------------------------------- @


if not _equations; @ ammend the linear policy rules if required and jump to the computation of impulse responses @
    if _ns>0;
        Lyx=Lyx|Lsx;
        Lyz=Lyz|Lsz;       
    endif;
    goto _DSGE_next1; 
endif;

_zstar=zeros(_nz,1);
_wstar=_xstar|_ystar|_zstar;
gmat=CDJac(&Sys,_wstar|_wstar,_nx+_ny);


if _GSchur;
    Amat=-gmat[.,1:_nx+_ny];
    Bmat= gmat[.,_nx+_ny+_nz+1:2*(_nx+_ny)+_nz];
    Cmat=-gmat[.,_nx+_ny+1:_nx+_ny+_nz]-gmat[.,2*(_nx+_ny)+_nz+1:2*(_nx+_ny+_nz)]*_Rho;
    {Lxx,Lxz,Llx,Llz}=SolveLA4(Amat,Bmat,Cmat,_Rho,_nx);
else;    
    _na=_nx+_ny+_nz;
    _nxy=_nx+_ny;
    Cu=gmat[1:_nu,1+_nx:_nx+_nu];
   Cxl=-gmat[1:_nu,1:_nx]~-gmat[1:_nu,1+_nx+_nu:_nxy];
    Cz=-gmat[1:_nu,1+_nx+_ny:_na];    
   Dxl=gmat[1+_nu:_nxy,1+_na:_na+_nx]~gmat[1+_nu:_nxy,1+_na+_nx+_nu:_na+_nxy];
   Fxl=gmat[1+_nu:_nxy,1:_nx]~gmat[1+_nu:_nxy,1+_nx+_nu:_nxy];
    Du=-gmat[1+_nu:_nxy,1+_na+_nx:_na+_nx+_nu];
    Fu=-gmat[1+_nu:_nxy,1+_nx:_nx+_nu];
    Dz=-gmat[1+_nu:_nxy,1+_na+_nxy:2*_na];
    Fz=-gmat[1+_nu:_nxy,1+_nxy:_na];
    {Lxx,Lxz,Llx,Llz,Lux,Luz}=SolveLA2(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,_Rho,_nx);      
    Lyx=Lux|Llx;
    Lyz=Luz|Llz;
endif;

output off;

if not _linear;
    @ Compute the Hessian Matrix @
    hcube=arrayinit((_nx+_ny)|(2*(_nx+_ny+_nz))|(2*(_nx+_ny+_nz)),0);
    for i (1,_nx+_ny,1);
        _eqno=i;
        hcube[i,.,.]=CDHesse(&Sys,_wstar|_wstar);
    endfor;

    {xcube,lcube}=SolveQA(Lxx,Lxz,Lyx,Lyz,gmat,hcube,_Rho,_Omega);

endif;

_DSGE_next1:

/* Compute Impulse Responses */
if _IR;
    struct PString {string name;};
    struct PString pstr;

    for i (1,_nz,1);
        if _equations; 
            {xt,lamt,zt}=Impulse2(nobs1,i);
        else;
            {xt,lamt,zt}=Impulse1(nobs1,i);
        endif;
        fig="Figure" $+ ftos(i,"*.*lf",1,0) $+ ".tkf";    
        Plott(fig,i);
    endfor;
endif;

/* Simulate the model and compute and output second moments */
if not _Moments; goto ende; endif;

if _equations;
    {sx,rx}=RBCRun2(nobs,nofs,_HPL);
else;
    {sx,rx}=RBCRun1(nobs,nofs,_HPL);
endif;

/* Output Results to the file outfile */

sml=GetMaxStringLength;
ccr=GetCrossCorrelations;
csx=GetRelSx;
if not ismiss(ccr); nccr=rows(ccr); else; nccr=0; endif;
if not ismiss(csx); ncsx=rows(csx); else; ncsx=0; endif;

cls;
outwidth 250;
output on;   
?"";
?"Second moments from simulated data:";
if _linear or _equations==0; ?"Linear policy functions were used"; else; ?"Quadratic policy functions were used"; endif;
?"";
for i (1,_nvar,1);
    if Var[i].varprint;
        vi=GetVarIndex(i);
        if _confidence;
            sx1=sortc(sx[vi,.]',1);
            i1=round(0.025*rows(sx1));
            i2=round(0.975*rows(sx1));
            str2=String1(sx1[i1],sx1[i2]);
            Ausstring2=strtruncpad(" ",sml) $+ str2;
            str1=String2(meanc(sx1));
            Ausstring1=strtruncpad(Var[i].varname,sml) $+ str1;
        else;
            Ausstring=strtruncpad(Var[i].varname,sml) $+ ftos(sx[vi],"*.*lf",8,2);
        endif;
        for j (1,ncsx,1);
            if _confidence;
                sx1=meanc(sx[vi,.]');
                sx2=meanc(sx[csx[j,2],.]');  
                str1=String2(sx1/sx2);
                Ausstring1=Ausstring1 $+ str1;
                Ausstring2=Ausstring2 $+ strtruncpad(" ",19);
            else;
                Ausstring=Ausstring $+ ftos(sx[vi]/sx[csx[j,2]],"*.*lf",8,2);
            endif;
        endfor;
        for j (1,nccr,1);
            if _confidence;
                nrx=cols(sx);
                rx1=zeros(nrx,1);                                          
                for i1 (1,nrx,1);
                    rx1[i1]=rx[i1,vi,ccr[j,2]];                    
                endfor;                               
                rx1=sortc(rx1,1);
                i1=round(0.025*nrx);
                i2=round(0.975*nrx);
                if vi==ccr[j,2];
                    str2=strtruncpad(" ",19);
                else;
                    str2=String1(rx1[i1],rx1[i2]);                
                endif;
                str1=String2(meanc(rx1));            
                Ausstring1=Ausstring1 $+  str1;
                Ausstring2=Ausstring2 $+  str2;                      
            else;
                Ausstring=Ausstring $+ ftos(rx[vi,ccr[j,2]],"*.*lf",8,2);
            endif;
         endfor; 
         if _confidence;
            for i1 (1,nrx,1);
                rx1[i1]=rx[i1,vi,_nvar+vi];
            endfor;
            rx1=sortc(rx1,1);
            str2=String1(rx1[i1],rx1[i2]);                
            str1=String2(meanc(rx1));            
            Ausstring1=Ausstring1 $+  str1;
            Ausstring2=Ausstring2 $+  str2;                
         else;       
            Ausstring=Ausstring $+ ftos(rx[vi,_nvar+vi],"*.*lf",8,2);
         endif;
         if _confidence;
            ?Ausstring1;
            ?Ausstring2;
         else;
            ?Ausstring;
         endif;         
    endif;
endfor;
?"";
?"Column 1: Variable name";
?"Column 2: Standard deviation";
for i (1,ncsx,1);
        Ausstring="Column " $+ ftos(2+i,"*.*lf",1,0) $+ ": Standard deviation relative to standard deviation of variable " $+ Var[csx[i,1]].varname;
        ?Ausstring;
endfor;
for i (1,nccr,1);
        Ausstring="Column " $+ ftos(2+ncsx+i,"*.*lf",1,0) $+ ": Cross correlation with " $+ Var[ccr[i,1]].varname;
        ?Ausstring;
endfor;
?"Column " ftos(3+ncsx+nccr,"*.*lf",1,0) ": First order autocorrelation";

output off;

/* Policy functions */

proc(1)=PF(typ,vnum,whut);

    local lhs;
    external matrix Lxx, Lxz, Lyx, Lyz, xcube, lcube, _linear;
    
    if typ=="x";
        lhs=(Lxx[vnum,.]~Lxz[vnum,.])*whut[1:rows(whut)-1];
        if _linear; retp(_xstar[vnum]+lhs); endif;
        lhs=lhs+0.5*(whut'*arraytomat(xcube[vnum,.,.])*whut);
        retp(_xstar[vnum]+lhs);
    elseif typ=="y";
        lhs=(Lyx[vnum,.]~Lyz[vnum,.])*whut[1:rows(whut)-1];
        if _linear; retp(_ystar[vnum]+lhs); endif;
        lhs=lhs+0.5*(whut'*arraytomat(lcube[vnum,.,.])*whut);
        retp(_ystar[vnum]+lhs);
    else;
        ?"Type must be either x or y"; ?"The program stops"; stop;
    endif;
    
endp;


/* This procedure uses the policy matrices Lxx, Lxz, Lyx, and Lyz and
** simulates the model and returns second moments.
*/

proc(2)=RBCRun1(nobs,nofs,hpl);

    local xt, zt, lamt, st, itn, epst, i, t,sx, rx, mom1,mom2, yt;
   
    if _confidence;
        sx=zeros(_nvar,nofs);
        i=nofs|2*_nvar|2*_nvar;
        rx=arrayinit(i,0);
    else;
        sx=zeros(_nvar,1);   
        rx=zeros(2*_nvar,2*_nvar);
    endif;

    if _loadShocks;
        {i,t}=fileinfo("eps_array.fmt");
        if i$=="";
            ?"Eps_array does not exist.";
            ?"Use _newShocks=1 or _MakesShocks=1 and restart the program";
            ?"The program stops now";
            stop;
        else;
            load eps_array;
        endif;
    endif;
      

  /* Simulation step */
  DosWin;
  cls; locate 5,5;?"Simulation Nr. ";
  for itn (1,nofs,1);
    locate 6,5;ftos(itn,"*.*lf",4,0);

    /* Initialize vectors */
  
    if _nx>0;
         xt=zeros(_nx,nobs+1);
    endif;
    lamt=zeros(_ny,nobs); 
    zt=zeros(_nz,nobs+1);
    if _loadShocks;
        epst=arraytomat(eps_array[itn,1:_nz,1:nobs+1]);
    else;
        epst=rndn(_nz,nobs+1);
    endif;  

    /* Compute time path */
    zt[.,1]=_Sigma*_Omega*epst[.,1];  
    for t (1,nobs,1);
         zt[.,t+1]=_Rho*zt[.,t]+_Sigma*_Omega*epst[.,t+1];         
         if _nx>0;
                xt[.,t+1]=Lxx*xt[.,t]+Lxz*zt[.,t];
                lamt[.,t]=Lyx*xt[.,t]+Lyz*zt[.,t];
         else;
                lamt[.,t]  =Lyx*xt[.,t]+Lyz*zt[.,t];
        endif;
    endfor;
    /* Filtering */
    st=zeros(nobs,_nvar);  
    for i (1,_nx,1);
          if HPL/=0; st[.,i]=HPFilter(xt[i,1:nobs]',HPL); else; st[.,i]=xt[i,1:nobs]';endif;
    endfor;
    for i (1,_ny,1);
          if HPL/=0; st[.,_nx+i]=HPFilter(lamt[i,.]',HPL); else; st[.,_nx+i]=lamt[i,.]';endif;
    endfor;
    for i (1,_nz,1);
        if HPL/=0;  st[.,_nx+_ny+i]=HPFilter(zt[i,1:nobs]',HPL); else; st[.,_nx+_ny+i]=zt[i,1:nobs]';endif;
    endfor;
  
    /* Computation of second moments */
    mom1=vcx(st);
    mom2=corrx(trimr(st,1,0)~trimr(st,0,1));
    if _confidence;
            sx[.,itn]=sqrt(diag(mom1));
            rx[itn,.,.]=mom2;
    else;
            sx=sx+sqrt(diag(mom1));
            rx=rx+mom2;
    endif;  
    
  endfor;
  DosWinCloseAll;
  if _confidence;
        retp(100*sx,rx);
  else;        
        retp((sx/nofs)*100,rx/nofs);
  endif;
endp;

/* This procedure uses the policy functions in PF from the main part
** simulates the model and returns second moments.
*/

proc(2)=RBCRun2(nobs,nofs,hpl);

    local xt, zt, lamt, st, itn, epst, i, t,sx, rx, mom1,mom2, yt;
   
    if _confidence;
        sx=zeros(_nvar,nofs);
        i=nofs|2*_nvar|2*_nvar;
        rx=arrayinit(i,0);
    else;
        sx=zeros(_nvar,1);   
        rx=zeros(2*_nvar,2*_nvar);
    endif;

    if _loadShocks;
        {i,t}=fileinfo("eps_array.fmt");
        if i$=="";
            ?"Eps_array does not exist.";
            ?"Use _newShocks=1 or _MakesShocks=1 and restart the program";
            ?"The program stops now";
            stop;
        else;
            load eps_array;
        endif;
    endif;
      

  /* Simulation step */
  DosWin;
  cls; locate 5,5;?"Simulation Nr. ";
  for itn (1,nofs,1);
    locate 6,5;ftos(itn,"*.*lf",4,0);

    /* Initialize vectors */
  
    if _nx>0;
         xt=zeros(_nx,nobs+1);
        xt[.,1]=_xstar;
    else;
        ?"Not implemented for _nx=0. The program stops";
        stop;
    endif;
    lamt=zeros(_ny,nobs); 
    zt=zeros(_nz,nobs+1);
    if _loadShocks;
        epst=arraytomat(eps_array[itn,1:_nz,1:nobs+1]);
    else;
        epst=rndn(_nz,nobs+1);
    endif;  

    /* Compute time path */
    zt[.,1]=_Sigma*_Omega*epst[.,1];  
    for t (1,nobs,1);
         zt[.,t+1]=_Rho*zt[.,t]+_Sigma*_Omega*epst[.,t+1]; 
        for i (1,_nx,1);
            xt[i,t+1]=PF("x",i,(xt[.,t]-_xstar|zt[.,t]|_sigma));
        endfor;
        for i (1,_ny,1);
           lamt[i,t]=PF("y",i,(xt[.,t]-_xstar|zt[.,t]|_sigma));
        endfor;
    endfor;

    /* Compute relative deviations */  
    for i (1,_nx,1);
        for t (1,nobs+1,1);
            xt[i,t]=(xt[i,t]-_xstar[i])/_xstar[i];
        endfor;
    endfor;
    for i (1,_ny,1);
        for t (1,nobs,1);
            lamt[i,t]=(lamt[i,t]-_ystar[i])/_ystar[i];
        endfor;
    endfor;
    /* Filtering */
    st=zeros(nobs,_nvar);  
    for i (1,_nx,1);
          if HPL/=0; st[.,i]=HPFilter(xt[i,1:nobs]',HPL); else; st[.,i]=xt[i,1:nobs]';endif;
    endfor;
    for i (1,_ny,1);
          if HPL/=0; st[.,_nx+i]=HPFilter(lamt[i,.]',HPL); else; st[.,_nx+i]=lamt[i,.]';endif;
    endfor;
    for i (1,_nz,1);
        if HPL/=0;  st[.,_nx+_ny+i]=HPFilter(zt[i,1:nobs]',HPL); else; st[.,_nx+_ny+i]=zt[i,1:nobs]';endif;
    endfor;
  
    /* Computation of second moments */
    mom1=vcx(st);
    mom2=corrx(trimr(st,1,0)~trimr(st,0,1));
    if _confidence;
            sx[.,itn]=sqrt(diag(mom1));
            rx[itn,.,.]=mom2;
    else;
            sx=sx+sqrt(diag(mom1));
            rx=rx+mom2;
    endif;  
    
  endfor;
  DosWinCloseAll;
  if _confidence;
        retp(100*sx,rx);
  else;        
        retp((sx/nofs)*100,rx/nofs);
  endif;
endp;

/* This procedure computes impulse reponse functions. It uses
** the linear policy functions in Lxx, Lxz, Lyx and Lyz.
**
** Impulses are to a one standard deviation shock in shock # i.
**
** The shock hit in period 2.
*/
proc(3)=Impulse1(nobs1,i);

    local xt,lamt,zt, t, evec;

    if _nx>0; xt=zeros(_nx,nobs1+1); else; xt=0; endif;    
  lamt=zeros(_ny,nobs1);
    zt=zeros(_nz,nobs1+1);
    evec=zeros(_nz,1);
    evec[i]=1;

    zt[.,2]=_Sigma*_Omega*evec;

    for t (2,nobs1,1);
        zt[.,t+1]=_Rho*zt[.,t];
        if _nx>0;
            xt[.,t+1]=Lxx*xt[.,t]+Lxz*zt[.,t];
           lamt[.,t] =Lyx*xt[.,t]+Lyz*zt[.,t];
        else;
          lamt[.,t]=Lyz*zt[.,t];
        endif;
    endfor;
    
    retp(real(xt),real(lamt),real(zt));
endp;


/* This procedure computes impulse reponse functions. It uses
** the policy functions in PF.
**
** Impulses are to a one standard deviation shock in shock # i.
**
** The shock hit in period 2.
*/
proc(3)=Impulse2(nobs1,i);

    local xt,lamt,zt, yt, t, evec;

    if _nx>0;
        xt=zeros(_nx,nobs1+1);
        xt[.,1]=_xstar;
        xt[.,2]=_xstar;
    else;
        ?"Not implemented for _nx=0. The program stops";
        stop;
     endif;    
    lamt=zeros(_ny,nobs1);
    lamt[.,1]=_ystar;
    lamt[.,2]=_ystar;
      zt=zeros(_nz,nobs1+1);    
     evec=zeros(_nz,1);
        evec[i]=1;

    zt[.,2]=_Sigma*_Omega*evec;

    for t (2,nobs1,1);
        zt[.,t+1]=_Rho*zt[.,t];
         for i (1,_nx,1);
             xt[i,t+1]=PF("x",i,(xt[.,t]-_xstar)|zt[.,t]|_sigma);
         endfor;
         for i (1,_ny,1);
            lamt[i,t]=PF("y",i,(xt[.,t]-_xstar)|zt[.,t]|_sigma);
         endfor;         
    endfor;
    @ Compute deviations @
    for i (1,_nx,1);
        for t (1,nobs1+1,1);
            xt[i,t]=(xt[i,t]-_xstar[i])./_xstar[i];
        endfor;
    endfor;
    for i (1,_ny,1);
        for t (1,nobs1,1);
            lamt[i,t]=(lamt[i,t]-_ystar[i])./_ystar[i];
        endfor;
    endfor;
    
    retp(real(xt),real(lamt),real(zt));
endp;

/* Plott: This procedure takes xt, lamt, yt, and zt from the
** main program and plots the impulse respones.
*/
proc(0)=Plott(fig,sn);

    local pmax, vmax, pmat, ord, ymin, ymax, ymaxall, yminall, i, j, k, t, pout;        

    {pmax,vmax}=GetMaxPlottNumber(sn);
    ord=(pmax|rows(lamt')|maxc(vmax));    
    ymax=zeros(pmax,1);
    ymin=zeros(pmax,1);
    k=zeros(pmax,1);    
    pstr=reshape(pstr, pmax, 1);    
    for i (1,pmax,1); pstr[i].name=" "; endfor; @ initialize legend strings @
    pmat=arrayinit(ord,0);

    for i (1,_nvar,1);     
        if Var[i].varplot;
            if Var[i].vartype=="z" and Var[i].varpos/=sn; continue; endif;  @ skip if not the current shock @
            j=var[i].plotno;            
            k[j]=k[j]+1;
            if Var[i].vartype=="x";
                pmat[j,.,k[j]]=xt[Var[i].varpos,1:cols(xt)-1]';
                pstr[j].name=pstr[j].name $+ Var[i].varname $+ "\000";                                        
            endif;
            if Var[i].vartype=="y";
                pmat[j,.,k[j]]=lamt[Var[i].varpos,.]';
                pstr[j].name=pstr[j].name $+ Var[i].varname $+ "\000";                                        
            endif;
            if Var[i].vartype=="z";                
                pmat[j,.,k[j]]=zt[Var[i].varpos,1:cols(zt)-1]';                
                pstr[j].name=pstr[j].name $+ Var[i].varname $+ "\000";
            endif;
        endif;
    endfor;    
    
    pmat=pmat*100;
    for i (1,pmax,1);
        ymax[i]=maxc(maxc(arraytomat(pmat[i,.,.])));
        ymin[i]=minc(minc(arraytomat(pmat[i,.,.])));
    endfor;
    ymaxall=maxc(maxc(ymax));
    yminall=minc(minc(ymin));
    GraphSettings;
    _ptek=fig;
    _plwidth=7;    
    
    if _nocolor; _pcolor={0,0,0,0,0,0,0}; _pltype={6,3,4,5,2,1}; endif;
    j=rows(lamt');
    if pmax==1;
        _plegstr=pstr[1].name;    
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        SetLegend(j,yminall,ymaxall,1,sn);
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");        
        pout=arraytomat(pmat[1,.,.]);
        xy(t,pout[.,1:vmax[1]]);
    endif;
    if pmax==2;
        begwind;
        window(2,1,0);
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        if _scale; SetLegend(j,yminall,ymaxall,1,sn); else; SetLegend(j,ymin[1],ymax[1],1,sn); endif;
        _plegstr=pstr[1].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[1,.,.]);        
        if _scale;  scale((0|t),(yminall|ymaxall)); endif;
        xy(t,pout[.,1:vmax[1]]);
        nextwind; 
        if _scale; SetLegend(j,yminall,ymaxall,2,sn); else; SetLegend(j,ymin[2],ymax[2],2,sn); endif;
        _plegstr=pstr[2].name;
        pout=arraytomat(pmat[2,.,.]);
        xy(t,pout[.,1:vmax[2]]);
        endwind;
    endif;
    if pmax==3;
        begwind;
        window(2,2,0);
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        if _scale; SetLegend(j,yminall,ymaxall,1,sn); else; SetLegend(j,ymin[1],ymax[1],1,sn); endif;
        _plegstr=pstr[1].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[1,.,.]);
        xy(t,pout[.,1:vmax[1]]);
        nextwind;        
        if _scale; SetLegend(j,yminall,ymaxall,2,sn); else; SetLegend(j,ymin[2],ymax[2],2,sn); endif;
        _plegstr=pstr[2].name;
        pout=arraytomat(pmat[2,.,.]);
        xy(t,pout[.,1:vmax[2]]);
        nextwind;
        _plegstr=pstr[3].name;
        if _scale; SetLegend(j,yminall,ymaxall,3,sn); else; SetLegend(j,ymin[3],ymax[3],3,sn); endif;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[3,.,.]);
        xy(t,pout[.,1:vmax[3]]);
        endwind;
    endif;    

   if pmax==4;        
        begwind;
        window(2,2,0);
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        if _scale; SetLegend(j,yminall,ymaxall,1,sn); else; SetLegend(j,ymin[1],ymax[1],1,sn); endif;
        _plegstr=pstr[1].name;          
        ?_plegstr;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[1,.,.]);
        if _scale;  scale((0|t),(ymin|ymax)); endif;
        xy(t,pout[.,1:vmax[1]]);
        nextwind;
        if _scale;SetLegend(j,yminall,ymaxall,2,sn); else; SetLegend(j,ymin[2],ymax[2],2,sn); endif;
        _plegstr=pstr[2].name;
        pout=arraytomat(pmat[2,.,.]);
        xy(t,pout[.,1:vmax[2]]);
        nextwind;
        if _scale;  SetLegend(j,yminall,ymaxall,3,sn); else; SetLegend(j,ymin[3],ymax[3],3,sn); endif;
        _plegstr=pstr[3].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[3,.,.]);
        xy(t,pout[.,1:vmax[3]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,4,sn); else; SetLegend(j,ymin[4],ymax[4],4,sn); endif;
        _plegstr=pstr[4].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[4,.,.]);
        xy(t,pout[.,1:vmax[4]]);
        endwind;
    endif;    

    if pmax==5;        
        begwind;
        window(3,2,0);
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        _pnumht=0.20;        
        _paxht=0.22;
        if _scale; SetLegend(j,yminall,ymaxall,1,sn); else; SetLegend(j,ymin[1],ymax[1],1,sn); endif;
        _plegstr=pstr[1].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[1,.,.]);
        if _scale;  scale((0|t),(ymin|ymax)); endif;
        xy(t,pout[.,1:vmax[1]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,2,sn); else; SetLegend(j,ymin[2],ymax[2],2,sn); endif;
        _plegstr=pstr[2].name;
        pout=arraytomat(pmat[2,.,.]);
        xy(t,pout[.,1:vmax[2]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,3,sn); else; SetLegend(j,ymin[3],ymax[3],3,sn); endif;
        _plegstr=pstr[3].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[3,.,.]);
        xy(t,pout[.,1:vmax[3]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,4,sn);else;SetLegend(j,ymin[4],ymax[4],4,sn); endif;
        _plegstr=pstr[4].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[4,.,.]);
        xy(t,pout[.,1:vmax[4]]);
        nextwind;
        if _scale;SetLegend(j,yminall,ymaxall,5,sn);else;SetLegend(j,ymin[5],ymax[5],5,sn); endif;
        _plegstr=pstr[5].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[5,.,.]);
        xy(t,pout[.,1:vmax[5]]);
        endwind;
    endif;    

    if pmax==6;
        begwind;
        window(3,2,0);
        t=seqa(1,1,j);
        _pline=(1~3~0~0~maxc(t)~0~1~0~7);
        _pnumht=0.20;        
        _paxht=0.22;
        if _scale; SetLegend(j,yminall,ymaxall,1,sn); else; SetLegend(j,ymin[1],ymax[1],1,sn); endif;
        _plegstr=pstr[1].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[1,.,.]);
        if _scale;  scale((0|t),(ymin|ymax)); endif;
        xy(t,pout[.,1:vmax[1]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,2,sn); else; SetLegend(j,ymin[2],ymax[2],2,sn); endif;
        _plegstr=pstr[2].name;
        pout=arraytomat(pmat[2,.,.]);
        xy(t,pout[.,1:vmax[2]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,3,sn); else; SetLegend(j,ymin[3],ymax[3],3,sn); endif;
        _plegstr=pstr[3].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[3,.,.]);
        xy(t,pout[.,1:vmax[3]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,4,sn); else; SetLegend(j,ymin[4],ymax[4],4,sn); endif;
        _plegstr=pstr[4].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[4,.,.]);
        xy(t,pout[.,1:vmax[4]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,5,sn); else; SetLegend(j,ymin[5],ymax[5],5,sn); endif;
        _plegstr=pstr[5].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[5,.,.]);
        xy(t,pout[.,1:vmax[5]]);
        nextwind;
        if _scale; SetLegend(j,yminall,ymaxall,6,sn); else; SetLegend(j,ymin[6],ymax[6],6,sn); endif;
        _plegstr=pstr[6].name;
        xlabel("Period");
        ylabel("Percentage Deviation from Trend");
        pout=arraytomat(pmat[6,.,.]);
        xy(t,pout[.,1:vmax[6]]);
        endwind;
    endif;    


    
endp;

@ ------------------------------------- MyDate -------------------------------

   Purpose: Returns the current date and time from
            Gauss intrinsic commands and prints it to an open file.
          
            Day Month Year
            Hours:Minutes:Seconds

          

   Usage: MyDate;

--------------------------------------------------------------------------- @
struct mstr {string m;};
struct mstr month;

proc(0)=MyDate;

  local t,d,str1,str2;
 
  month=reshape(month,12,1);

  month[1].m="January";
  month[2].m="Feburary";
  month[3].m="March";
  month[4].m="April";
  month[5].m="May";
  month[6].m="June";
  month[7].m="July";
  month[8].m="August";
  month[9].m="September";
  month[10].m="October";
  month[11].m="November";
  month[12].m="December";

  d=date;
  t=time;
  str1=ftos(d[3],"%0*.*lf",2,0) $+ " " $+ month[d[2]].m $+" " $+ ftos(d[1],"%0*.*lf",2,0);
  str2=ftos(t[1],"%*.*lf",2,0)  $+ ":" $+ ftos(t[2],"%0*.*lf",2,0) $+":" $+ ftos(t[3],"%0*.*lf",2,0);
  ?str1;
  ?str2;
retp;

endp;

@ ------------------------------- HPFilter ---------------------------------------------

   Purpose:  returns the cyclical component of a time series x
             as defined by the Hodrick-Prescott filter.

   usage: c=HPFilter(x,mu);

   Input:   x : nobs x 1, vector with observations

           mu : scalar, the filter weight (i.e. mu=1600 for quarterly data)

   Output   c : nobs x 1, vector, the cyclical component of x
------------------------------------------------------------------------------------------ @


proc(1)=HPFilter(x,mu);

local m, i, nobs;

  nobs=rows(x);
  
  m=zeros(nobs,3);

  /* Set up the matrix */
  i=1;
  do until i>nobs;
   m[i,1]=mu;
   m[i,2]=-4*mu;
   m[i,3]=1+6*mu;
        i=i+1;
  endo;
  m[1,1]=0;         m[2,1]=0;
  m[1,2]=0;         m[2,2]=-2*mu;     m[nobs,2]=m[2,2];
  m[1,3]=1+mu;      m[2,3]=1+5*mu;
  m[nobs,3]=m[1,3]; m[nobs-1,3]=m[2,3];
 
 retp(x-bandsolpd(x,m));  @ The command bandlospd provides the solution of equation (A.4.2) @

endp;   


@ ------------------------------- GraphSettings ---------------------------------------

   Purpose: Changes Gauss' default initialization of graphic routines

   Usage: GraphSettings;
---------------------------------------------------------------------------------------- @

proc(0)=GraphSettings;

 external matrix  _pltype, _ptitlht, _pnumht,
         _paxht, _pmcolor, _plwidth, _pcsel, _pdate;
 external proc GraphSet;
 GraphSet;
   _pdate="";        @ do not plot date               @
  _pltype=6;         @ solid lines only               @
 _plwidth=4;         @ line width always 4            @
   _paxht=0.15;      @ size of axis labels            @        
  _pnumht=0.15;      @ size of axis numbering         @
 _ptitlht=0.20;      @ size of titles                 @
 _pmcolor={0,        @ color of axis: black           @
           0,        @ color of axis' numbers         @
           0,        @ color of x-axis label          @
           0,        @ color of y-axis label          @
           0,        @ color of z-axis label          @
           0,        @ color of title                 @
           0,        @ color of boxes                 @
           0,        @ color of date and author       @
           15};      @ background color:        white @
   _pcsel={0,        @ color of first line: black     @
           9,        @ color of second line: light blue @
           10,       @ color of third line:  light green @
           12,       @ color of fourth line: light red      @
           13,       @ color of fifth line:  light magenta @
           11,       @ color of sixth line: light cyan @
           6};       @ color of seventh line: brown @
   fonts("Simplex, Simgrma");
   _plctrl={0,0,0,0,0,0,1,1,1}; @ lines 7-9 with symbols @
   _pstype={1,1,1,1,1,1,8,9,10}; @ lines 7-9: solid circle, solid square, solid triangle @
   _psymsiz=2;          
    retp;
endp;           

@ ------------------------------------ CDJac --------------------------------


  purpose: computes a central difference approximation of the Jacobian
          matrix of a system of n non-linear functions y_i=f^i(x), where
          x is a column vector of dimension m.

  usage: Jac=CDJac(&f,x0,n)


 
 i nput:  &f: pointer to the routine that returns the m-vector f(x)
         x0: m vector x0, the point at which the derivatives are to be evaluated
         n : the size of the vector f(x)

  output: Jac: n by m matrix of parital derivatives

  algorithm: based on (A.2.8) in Heer and Maussner, see also Dennis and Schnabel (1983), 
            Algorithm A5.6.4.
--------------------------------------------------------------------------------------- @

proc(1)=CDJac(&f,x0,n);

  local eps, i, j, m, h, df, f1, f2, x1, x2, temp;
  local f:proc;
  
  m=rows(x0);
  df=zeros(n,m);
  eps=MachEps^(1/3);
  x1=x0;
  x2=x0;
 
  for i (1, m, 1);
     if x0[i]<0;
        h=-eps*maxc(abs(x0[i])|1.0);
     else;
        h=eps*maxc(abs(x0[i])|1.0);
    endif;
    temp=x0[i];
    x1[i]=temp+h;
    x2[i]=temp-h;
        h=x1[i]-temp; @ Trick to increase precision slightly, see Dennis and Schnabel (1983), p. 99 @
    f1=f(x1);
    f2=f(x2);
    for j (1,n,1);    
       df[j,i]=(f1[j]-f2[j])/(2*h);       
    endfor;
    x1[i]=x0[i];
    x2[i]=x0[i];    
 endfor;
 retp(df);

endp;

/* CDHesse
**
** usage Hmat=CDHesse(&f,x0)
**
** purpose: compute approximation of the Hesse matrix of a real valued function
**          y=f(x), where x is a n vector
**
** input:  &f: pointer to the procedure that returns f(x)
**         x0: n vector, the point at which the Hesse matrix is to be computed
**
** output: Hmat: n by n matrix
**
** algorithm: based on the formulas (A.2.10) and (A.2.11) in Heer and Maussner
*/


proc(1)=CDHesse(&f,x0);

 local i, j, n, hmat, h, eps, temp, x1, x2,x3,x4;
 local f:proc;

 n=rows(x0);
 h=zeros(n,1);
 hmat=zeros(n,n);
 eps=MachEps^(1/3);
 for i (1,n,1);
    if x0[i] <0;  h[i]=-maxc(abs(x0[i])|1.)*eps; else; h[i]=maxc(abs(x0[i])|1.0)*eps; endif;
    temp=x0[i]+h[i];
    h[i]=temp-x0[i]; @ the last two steps increase precision slightly, see Dennis and Schnabel (1983), p. 321 for this trick @    
 endfor;

 for i (1,n,1);
     x1=x0;
     x2=x0;
  x1[i]=x1[i]+h[i];
  x2[i]=x2[i]-h[i];
    hmat[i,i]=(f(x1)+f(x2)-2*f(x0))/((h[i]^2));
    for j (i+1,n,1);    
           x1=x0; x1[i]=x1[i]+h[i]; x1[j]=x1[j]+h[j];
           x2=x0; x2[i]=x2[i]-h[i]; x2[j]=x2[j]+h[j];
           x3=x0; x3[i]=x3[i]+h[i]; x3[j]=x3[j]-h[j];
           x4=x0; x4[i]=x4[i]-h[i]; x4[j]=x4[j]-h[j];
           hmat[i,j]=(f(x1)-f(x2)-f(x3)+f(x4))/(4*h[i]*h[j]);
           hmat[j,i]=hmat[i,j];           
    endfor;
 endfor;

retp(hmat);

endp;


@ ------------------------- MachEps ------------------------------

  Purpose: Computes the machine epsilon
  Usage: e=MachEps;

  Output: the smallest number e so that (1+e)>1 is true

------------------------------------------------------------------- @

proc(1)=MachEps;

    local eps, i;
        eps=1;
        do while (1+eps)>1;
            eps=eps/2;
        endo;
    retp(2*eps);
endp;



/* The remaining procedures are used to output data. They are
** of no special interest and therefore I do not comment them.
*/

proc(2)=GetMaxPlottNumber(sn);
    local n0, n1, n2;
    n0=zeros(6,1);
    n1=0;
    for i (1,_nvar,1);
        if Var[i].varplot;
            if Var[i].vartype=="z" and Var[i].varpos/=sn; continue; endif;
            n0[var[i].plotno]=n0[var[i].plotno]+1;
            n2=var[i].plotno;
            if n2>n1;
                n1=n2;
            endif;
        endif;
    endfor;
    n2=maxc(n0);    
    retp(n1,n0);
endp;


proc(1)=GetMaxStringLength;

    local s1, s2, i;

    s1=strlen(Var[1].varname);
    for i (2,_nvar,1);
        s2=strlen(Var[i].varname);
            if s2>s1; s1=s2; endif;
    endfor;

    retp(s1);

endp;

proc(1)=GetCrossCorrelations;

    local i, m;
    m=(0~0);
    for i (1,_nvar,1);
        if Var[i].crosscorr;             
            if Var[i].vartype=="x"; m=m|(i~(Var[i].varpos)); endif;
            if Var[i].vartype=="y"; m=m|(i~(_nx+Var[i].varpos)); endif;
            if Var[i].vartype=="z"; m=m|(i~(_nx+_ny+Var[i].varpos)); endif;
        endif;
    endfor;
    if rows(m)>1;
        retp(m[2:rows(m),.]);
    else;   
        retp(miss(1,1)~miss(1,1));
    endif;

endp;

proc(1)=GetRelSx;

    local i, m;
    m=(0~0);
    for i (1,_nvar,1);
        if Var[i].relsx;             
            if Var[i].vartype=="x"; m=m|(i~(Var[i].varpos)); endif;
            if Var[i].vartype=="y"; m=m|(i~(_nx+Var[i].varpos)); endif;            
            if Var[i].vartype=="z"; m=m|(i~(_nx+_ny+Var[i].varpos)); endif;
        endif;
    endfor;
    if rows(m)>1;
        retp(m[2:rows(m),.]);
    else;
        retp(miss(1,1)~miss(1,1));
    endif;

endp;

            
proc(1)=GetVarIndex(i);

    local vi;

    if Var[i].vartype=="x"; vi=Var[i].varpos;         endif;
    if Var[i].vartype=="y"; vi=_nx+Var[i].varpos;     endif;    
    if Var[i].vartype=="z"; vi=_nx+_ny+Var[i].varpos; endif;

    retp(vi);

endp;


proc(0)=SetLegend(j,ymin,ymax,i,sn);
    local a,b;
    external matrix _plegctl, _plegend;

    if ismiss(_plegend[sn,i,1]);
        if abs(ymin)>ymax;
           _plegctl=(1~6~floor(j/3)~ymin);
        else;
          _plegctl=(1~6~floor(j/3)~0.5*ymax);
        endif; 
    else;     
       a=arraytomat(_plegend[sn,i,1]);
       b=arraytomat(_plegend[sn,i,2]);
        _plegctl=(1~6~a~b);
       
    endif;      
    
retp;

endp;

proc(1)=String1(s1,s2);

    local c;    
    c="  (" $+ ftos(s1,"*.*lf",6,2) $+ "-" $+ ftos(s2,"*.*lf",6,2) $+ ")  ";
    
retp(c);

endp;

proc(1)=String2(s);
    local l1, i1, s1;
    s1=ftos(s,"*.*lf",6,2);
    l1=strlen(s1);
    i1=round((19-l1)/2);
    s1=strtruncpad(" ",i1) $+ s1;
    i1=strlen(s1);
    s1=s1 $+ strtruncpad(" ",19-i1);
 retp(s1);
endp;


#include CoRRAM_4.src;
#include CoRRAM_5.src;
