/*               vsimul.prg                                               */
/* simulates two-component ar(1) model for surplus, v
   run rettos.prg first                                                   */
/* ---------------------------------------------------------------------- */

output file = vsimul.out reset;

library pgraph;
bandw = 1;
convt = 0;
graphs = 1;

/* --------- */
/* LOAD DATA */
/* --------- */

/* uses s, rgov, matstr2, price, money base, fama bliss
   loaded from rtos.prg  This is straight from cambpl2.prg
*/

rs = s./apunew[2:rows(apunew)]*100;
                                /* NOTE price level is relative to 1, not 100.
                                   you MUST keep this convention below, use
                                   0.5 not 50 for price level */


T = rows(rs);                       /* rows of changes -- 61-96*/
Tl = T+1;                           /* rows of levels -- 60-96 */
dates = seqa(60,1,Tl);              /* NOTE 123150 is graphed as 50, not 51 */

/* ------------------------------ */
/* now load GDP, consumption data */
/* ------------------------------ */

load x[] = gdp.dat;
x = reshape(x,rows(x)/4,4);
yrq2 = x[.,1];
qq2 = x[.,2];
gdpq = x[.,4];

load x[] = cons.dat;
x = reshape(x,rows(x)/10,10);
yr3   = x[.,1];
qq3   = x[.,2];
GCNQ = x[.,7];
GCSQ = x[.,10];

gdpa = chgfreq(gdpq,1,4,0);    /* option to use 4th quarter values */
ca = chgfreq(gcnq+gcsq,1,4,0);

gdpa= gdpa[2:rows(gdpa)]; /* levels so */
ca = ca[2:rows(ca)];      /* chop to 60:96 */

/* ------------------------------------ */
/* plot variables to check stationarity */
/* ------------------------------------ */

rmkt2t = (mkt2+abase)./apunew*100;       /* real market value incl. m base */
nmkt2t = mkt2+abase;
stoct = ((rs/1000)./ca[2:Tl]);
rrgovt = rgov./infl;
vtoct = ((rmkt2t/1000)./ca);
vntoct = (nmkt2t./1000)./ca;
dct   = (ca[2:Tl]./ca[1:Tl-1]);
/* graph basic series */

if graphs;

    dv = dif(ln(rmkt2t));
    dvn = dif(ln(nmkt2t));
    "Real debt growth on nominal debt growth";
    __altnam = 0; call ols("",dv,dvn);
    graphjc;
    title("real and nominal debt growth");
    xy(dates[2:rows(dates)],dv~dvn)



    graphjc;
    title("Debt, Surplus and Inflation");
    if bandw;
        _pcolor = 15;
    endif;
    _plegstr = "v/c - 45\000s/c\000inflation";
    _plegctl = 1;
    _pltype = 3|6|6;
    _plctrl = 0|1|0;
    _psymsiz = 2;
     xtics(60,96,4,4);
    ytics(-16,16,4,4);
    px = fill(dates,dates[2:rows(dates)].*(1~1));
    py = fill(100*(vtoct-0.45),
              100*(stoct)~
              100*(ln(infl)) );

        /*       1 typ xs ys xe ye 1=new co thick(0) */
        _pline = 1~6~61~0~97~0~1~15~0;
    xy(px,py);
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim1.pic;
    endif;


endif;

/* ------------ */
/*  Now run VAR */
/* ------------ */

/* Simple VAR */

x = stoct~vtoct[2:rows(vtoct)];
x = x-meanc(x)';
T = rows(x);
A = x[2:T,.]/x[1:T-1,.];
e = x[2:T,.]-x[1:T-1,.]*A;
vmat = e'e/T;
corrm = vmat./((diag(vmat)*diag(vmat)')^0.5);
A = A';
"";
"Simple var"   ;
"Coefficients of sc vc on one lag";
"  sc(t) " A[1,.];
"  vc(t) " A[2,.];
"";
"correlation matrix";
corrm;
"standard deviations";
diag(vmat)^0.5;

"";

"OLS version of simple VAR";
let __altnam = "const" "sc(t-1)" "vc(t-1)" "sc(t)" ;
call ols("",x[2:T,1],x[1:T-1,.]);

let __altnam = "const" "sc(t-1)" "vc(t-1)" "vc(t)" ;
call ols("",x[2:T,2],x[1:T-1,.]);

"variance /covariance matrix stored in vmat";

/* inflation in the simple VAR. DO we need another state variable?*/

"inflation on contemporaneous these variables";
y = ( ln(infl) -meanc(ln(infl)) );
alphols = ( y/x)';
"  pi(t) ";; alphols;
"OLS version ";
let __altnam = "const" "sc(t)" "vc(t)" "pi(t)" ;
call ols("",y,x);

"subsample coefficients";
"full sample " (y/x)';
"      60-70 " (y[1:10]/(x[1:10,.]~ones(10,1)))';
"      70-83 " (y[11:23]/(x[11:23,.]~ones(13,1)))';
"      83-96 " (y[23:36]/(x[23:36,.]~ones(14,1)))';

/*
/* Simple VAR with inflation */

x = x~( ln(infl) -meanc(ln(infl)) );
T = rows(x);
A1 = x[2:T,.]/x[1:T-1,.];
e1 = x[2:T,.]-x[1:T-1,.]*A1;
vmat2 = e1'e1/T;
corrm2 = vmat2./((diag(vmat2)*diag(vmat2)')^0.5);
A1 = A1';
"";
"Simple var"   ;
"Coefficients of vc sc on one lag";
"  sc(t) " A1[1,.];
"  vc(t) " A1[2,.];
"infl(t) " A1[3,.];
"";
"correlation matrix";
corrm2;

"variance /covariance matrix stored in vmat2";

"OLS version of simple VAR with inflation";

let __altnam = "const" "sc(t-1)" "vc(t-1)" "pi(t-1)" "sc(t)";
call ols("",x[2:T,1],x[1:T-1,.]);

let __altnam = "const" "sc(t-1)" "vc(t-1)" "pi(t-1)" "vc(t)";
call ols("",x[2:T,2],x[1:T-1,.]);

let __altnam = "const" "sc(t-1)" "vc(t-1)" "pi(t-1)" "pi(t)";
call ols("",x[2:T,3],x[1:T-1,.]);

*/


/* ---------------------- */
/* SIMULATION             */
/* ---------------------- */

/* input A in s,v process. only a11 is used to find b, a22 must = a11 */


/* Use OLS A, but impose identifying restriction */
A[2,1] = A[1,1];

/* Direct:
A = ( 0.5 ~ 0.06 )|
    (-0.5 ~ 0.96 );
*/
/* input how inflation depends on s,v */
        /* s     v   */
alphsv = (1 ~ -0.216947 );

/* input V in s,v, process via sigma and rho */

/* Option to use OLS residuals */
/*
sigs = diag(vmat)^0.5;
sigv = sigs[2];
sigs = sigs[1];
rhosv = vmat[2,1]/(sigs*sigv);
*/
/* create residuals from s,v ar1 using restricted model */

x = stoct~vtoct[2:rows(vtoct)];
x = x-meanc(x)';
svresid = x[2:rows(x),.] - x[1:rows(x)-1,.]*A';
svmat = svresid'svresid/rows(svresid);
"10000* restricted residual vcov matrix" 10000*svmat;
"Stand dev's, correlation" svmat[1,1]^0.5 svmat[2,2]^0.5
                            svmat[1,2]/(svmat[1,1]*svmat[2,2])^0.5;
olsmat = e'e/rows(e);
"10000* OLS residual vcov matris       " 10000*olsmat;
"std devs and correlation"  olsmat[1,1]^0.5 olsmat[2,2]^0.5
                            olsmat[1,2]/(olsmat[1,1]*olsmat[2,2])^0.5;

/* use residuals from restricted model for error term*/
sigs = svmat[1,1]^0.5;
sigv = svmat[2,2]^0.5;
rhosv = svmat[1,2]/(svmat[1,1]*svmat[2,2])^0.5;


proc(7) = dosim(A,alphsv,sigv,sigs,rhosv,svresid);
    /* if svresid == T, will simulate from random normal T observations */
    /* test is on rows of svresid; if rows > 1 treats as residual series */


    /* also uses data from above passed as globals! */

    local eta,gam,bet,B,D,V,cv,DV,sigmat,sigz,sigc,rhocz,C,Emat,alphcz,vcv,svt1,
    scsim,vcsim,dpsim,pusim,vncsim,dvncsim,dvcsim,j,maxit;

    /* find coefficients in c,z process of s components */

    eta = (a[2,2]+a[1,1] + ((a[2,2]-a[1,1])^2 - 4*a[1,2]*a[1,1])^0.5 )/2;
    gam = (a[2,2]+a[1,1] - ((a[2,2]-a[1,1])^2 - 4*a[1,2]*a[1,1])^0.5 )/2;
    bet = 1/(a[2,2]+a[1,2]);

    B = ( 1 ~ 1 )|
        ( (bet*eta/(1-bet*eta)) ~ (bet*gam/(1-bet*gam)) );

    D = (eta~0)|(0~gam);

    alphcz = alphsv*inv(B);

    V = ( ( sigs^2 )~(sigv*sigs*rhosv) )|
        ( (sigv*sigs*rhosv) ~ (sigv^2) );

    CV = chol(v)';

    sigmat = inv(B)*V*inv(B');
    sigz = sigmat[1,1]^0.5;
    sigc = sigmat[2,2]^0.5;
    rhocz = sigmat[1,2]/(sigmat[1,1]*sigmat[2,2])^0.5;

    C = chol(sigmat)'; /* now cc' = sigmat ; c*normal(0,1) is error*/

    /* print */
    "";
    "Calling Simulation";
    "";
    "Sssumed s,v process ";
    "A s,v on lagged s,v "; A;
    "sigs sigv rho       " sigs sigv rhosv;
    "Dp on current s,v   " alphsv;

    "";
    "Implied c,z process ";
    "diagonals of D " eta gam;
    "beta           " bet;
    "sigc sigz rho  " sigc sigz rhocz;
    "Dp on curr c,z " alphcz;

    /* give interpretation in terms of  v,s-v process */

    Emat =  (1~(-(1-bet)/bet))|
             (0~1);

    "";
    "Implied s-rv, v process ";
    "AR(1) matrix" Emat*A*inv(Emat);
    "sigma s-rv, sigma v and correl";
    vcv = Emat*V*Emat'; (diag(vcv)^0.5)' vcv[2,1]/(vcv[1,1]*vcv[2,2])^0.5;
    "loading of Dp on s-rv,v";

    alphsv*inv(Emat);


    /* form time series*/

    if rows(svresid) > 1;
        "Using residuals";
         maxit = rows(svresid);
    else;
        "generating artificial time series";
        maxit = svresid; maxit;
    endif;

    j = 1;
    svt1 = (stoct[1]-meanc(stoct))~(vtoct[2]-meanc(vtoct[2:rows(vtoct)]));

    do while j <= maxit;
        if rows(svresid) >  1;     /* use residuals from above */
            svt1 = svt1 | (( (A*svt1[j,.]') + svresid[j,.]' )');
        else;
            svt1 = svt1 | (( (A*svt1[j,.]') + (-cv*rndn(2,1)) )');
        endif;
        j = j+1;
    endo;

    scsim = svt1[.,1]+meanc(stoct);
    vcsim = svt1[.,2]+meanc(vtoct[2:rows(vtoct)]);
    dpsim = svt1*alphsv'+meanc(ln(infl));

    /* apunew = apunew[1]*(1|cumprodc(infl));*/

    pusim  = apunew[1]*(1|exp(cumsumc(dpsim)));
    vncsim = vntoct[1]|(vcsim.*(pusim[2:rows(pusim)]./100));
    dvncsim = dif(ln(vncsim));
    dvcsim = ln(vtoct[2]/vtoct[1])|dif(ln(vcsim));

    retp(scsim,vcsim,dpsim,pusim,vncsim,dvncsim,dvcsim);
endp;


/* do simulation with artificial data */

rndseed 2010;
{scsim,vcsim,dpsim,pusim,vnsim,dvncsim,dvcsim} =
            dosim(A,alphsv,sigv,sigs,rhosv,rows(svresid));

if graphs;
    graphjc;

    graphjc;
    fonts("simplex simgrma");
    title("Simulated surplus and debt growth");
    if bandw;
        _pcolor = 15;
    endif;
    _plegstr = "surplus/value\000value growth";
    _plegctl = 1;
    _pltype = 3|6|6;
    _plctrl = 0|1|0;
    _psymsiz = 2;
     xtics(60,96,4,4);
    ytics(-16,16,4,4);
        /*       1 typ xs ys xe ye 1=new co thick(0) */
        _pline = 1~6~61~0~97~0~1~15~0;
    xy(dates[3:rows(dates)],
      (100*(vcsim[2:rows(vcsim)]./vcsim[1:rows(vcsim)-1]-1))~
      (100*(scsim[2:rows(scsim)]./vcsim[1:rows(vcsim)-1])));
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim2a.pic;
    endif;


    fonts("simplex simgrma");
    title("Simulated v/c, s/c and Inflation");
    if bandw;
        _pcolor = 15;
    endif;
    _plegstr = "v/c - 45\000s/c\000inflation \202D\201p";
    _plegctl = 1;
    _pltype = 3|6|6;
    _plctrl = 0|1|0;
    _psymsiz = 2;
     xtics(60,96,4,4);
    ytics(-16,16,4,4);
        /*       1 typ xs ys xe ye 1=new co thick(0) */
        _pline = 1~6~61~0~97~0~1~15~0;
    xy(dates[2:rows(dates)],100*((vcsim-0.45)~scsim~dpsim));
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim2.pic;
    endif;
endif;



/* do simulation with residuals */

{scsim,vcsim,dpsim,pusim,vnsim,dvncsim,dvcsim} =
            dosim(A,alphsv,sigv,sigs,rhosv,svresid);

/* see if OLS gives closer inflation prediction */

{scsim3,vcsim3,dpsim3,pusim3,vnsim3,dvncsim3,dvcsim3} =
            dosim(A,alphols,
@(0.1~-0.216947),@
sigv,sigs,rhosv,svresid);

if graphs;
    graphjc;
    if bandw;
        _pcolor = 15;
    endif;
    title("Actual and Simulated Inflation\L"$+
                        "\202D\201p = 1.0 sc - 0.21 vc") ;
    fonts("simplex simgrma");
    _plegstr = "Actual\000Simulated" ;
    _plegctl = 1|4|84|8;
    _psymsiz = 2;
    _pltype = 6|3;
    _plctrl = 0|0;
    ytics(-2,14,2,2);
    ylabel("Percent inflation");
    /*       1 typ xs ys xe ye 1=new co thick(0) */
    _pline = 1~6~61~0~97~0~1~15~0;
    xy(dates[2:rows(dates)],100*(ln(infl)~dpsim));
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim3a.pic;
    endif;

    title("Actual and Simulated Inflation\L"$+
                        "\202D\201p = 0.08 sc - 0.21 vc (OLS)") ;
    xy(dates[2:rows(dates)],100*(ln(infl)~dpsim3));
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim3b.pic;
    endif;
endif;

"std of actual - simulated inflation ";
100*stdc(ln(infl)-dpsim);
"std of actual - simulated inflation, OLS d ";
100*stdc(ln(infl)-dpsim3);


/* study real and nominal debt growth, and zero inflation debt growth */
if graphs;
    graphjc;
    if bandw;
        _pcolor = 15;
    endif;
    _plegstr = "Real debt growth"$+
              "\000Simulated nominal debt growth\000Actual nominal debt growth";
    _plegctl = 1;
    _pltype = 6|6|3;
    _plctrl = 0|1|0;
    _psymsiz = 2;
     /*       1 typ xs ys xe ye 1=new co thick(0) */
     _pline = 1~6~61~0~97~0~1~15~0;
    xtics(60,96,4,4);
    ytics(-12,20,4,4);
    ylabel("Percent growth");
    xy(dates[2:rows(dates)],
    100*(@dvcsim~@dif(ln(vtoct))~dvncsim~dif(ln(vntoct))));
    if convt;
        dos pqgrun graphic.tkf -c=2 -cf=sim4.pic;
    endif;
endif;






output off;