/* ------------------------------------------------------------------------- */
/*                            pihsim.gpg                                     */
/*  simulates pih model, runs fama regressions                               */
/* ------------------------------------------------------------------------- */

library pgraph;
checkanl = 0;                                 /* check analytical formulas? */
dographs = 1;                                 /* make graphs of simulation  */
graphcig = 1;                                 /* make graph of level of c,i,g*/
graphrat = 1;                                 /* graphs of cointeg., y-c, etc*/
output file=pihsim.out reset;

r     = .1;
delta = .1;
rho   = .8;
sige  = 5;                    /* for ebar=100, pick % std dev of labor income */
sigeps = (1-rho^2)^(1/2);     /* std(e) = 1/(1-rho^2)^(1/2)std(eps)           */
ebar  = 100;
capT  = 100;

theta = r+delta;                                        /* r = theta - delta */
k1 = ebar/(2*theta);                    /* certainty: y = theta*k + e ,      */
                                        /* want theta*k = 33%, e = 66%       */
                                        /* thus theta k = e/2 , k = e/2*theta*/
lam  = 1/(1+r);
m    = r/(1+r-rho);
gam  = (1+rho*r-rho)/(1+r-rho) ;
eta  = (1+rho*r-rho-delta*(1-rho))/(1+r-rho);
sqig = (1-rho-delta)/(1+r-rho);

/* draw e */
eps  = sige*rndn(capT,1);
etil =  recserar(eps,0,rho);
e = ebar + etil;

/* find consumption and capital */
c = zeros(capT,1);
k = zeros(capT,1);
k[1] = k1;

t = 1;
do while t <= capT;
    c[t] = r*k[t] + ebar + m*etil[t];
    if (t /= capT);
        k[t+1] = (1+r)*k[t] + e[t] - c[t];
    endif;
    t = t+1;
endo;

/* find other series */

yn = r*k + e;
y  = (r+delta)*k + e;
ni = yn - c;
gi = ni + delta*k;
w  = c./r;

/* check analytical formulas */
if checkanl;
     "roots eta and gamma" eta gam;
     "max divergence between analytical and constructed ";
     "DIfferences:" ;
    dc  = c[2:capT]-c[1:capT-1];
    "    consumption" maxc(dc - m*eps[2:capT]);
    dw  = w[2:capT]-w[1:capT-1];
    "         Wealth" maxc(dw - m/r*eps[2:capT]) ;
    dk  = k[2:capT]-k[1:capT-1];
    "        Capital" maxc(dk - (1-m)*etil[1:capT-1]) ;
    dgi = gi[2:capT]-gi[1:capT-1];
    "     Gross inv."
                maxc(dgi - (1-m)*(etil[2:capT] - (1-delta)*etil[1:capT-1])) ;
    dyn = yn[2:capT]-yn[1:capT-1];
    "     Net income" maxc(dyn - (etil[2:capT] - gam*etil[1:capT-1])) ;
    dy  = y[2:capT]-y[1:capT-1];
    "            GNP" maxc(dy - (etil[2:capT] - eta*etil[1:capT-1])) ;
    "Levels: ";
    "         Wealth" maxc(w  - c./r);
    "        Capital" maxc(k  - (c./r - ebar/r - m/r*etil));
    "     Gross inv." maxc(gi - (delta/r*c - delta/r*ebar + sqig*etil));
    "            NNP" maxc(yn - (c+(1-m)*etil));
    "            GNP" maxc(y - ((1+delta/r)*c - delta/r*ebar + sqig*etil));
endif;


if dographs;
    graphset;
    _Pltype = 6;
    _Pframe = 0|0;
    _Pdate = "";
    _Ptek = "" ;
    tim = seqa(1,1,capT);
    if graphcig;
        xy(tim,(c~y~gi));
    endif;

    if graphrat;
         _Ptitlht = .2;
         xy(tim, (y - ((1+delta/r)*c - delta/r*ebar)));
    endif;


endif;

    format /RD 10,3;
    "  r " r    " delta" delta  "theta" theta;
    " rho" rho  "s(eps)" sigeps " s(e)"  sige;
    "ebar" ebar "    k0" k1     "    m"  m   ;
    " eta" eta  " gamma" gam    " sqig" sqig ;
    " yt = a(L)et, a(j) and partial sums of a(j) ";
     a = rho-eta;
     a = rho^(seqa(0,1,20))*a;
     a = 1|a;

     b = rho-(1-delta);
     b = rho^(seqa(0,1,20))*b;
     b = 1|b;

    let label = " j" "GNP a(j)" "part sum" "GI a(j)" "part sum";
    $label'
     seqa(0,1,21)~a~recserar(a,a[1],1)~b~recserar(b,b[1],1);

output off;