/* runvar.prg runs vars for ep article */

library pgraph;
convt = 0;
output file = runvar.out reset;

maxsim = 24;      /* length of irs */
k = 4;           /* number of lags in regressions */
levels = 0;      /* run in log levels or specify diff, ci structure */

load x[] = marty.dat;
x = reshape(x,rows(x)/8,8);
y = x[.,1];
p = x[.,2];
pcom = x[.,3];
nbr = x[.,4];
tr = x[.,5];
ff = x[.,6];
m1 = x[.,7];
m2 = x[.,8];
dates = seqb(65.5,95.25,0.25);

load x2[] = citirate.dat;
x2 = reshape(x2,rows(x2)/5,5);
cyr = x2[.,1];
cmo = x2[.,2];
fygl = x2[.,3];
fygn3 = x2[.,5];
fygl = chgfreq(fygl,1,3,0);     /* use last obs in qtr for qtr data */
fygn3 = chgfreq(fygn3,1,3,0);
cyr = chgfreq(cyr,1,3,0);
cmo = chgfreq(cmo,1,3,0);
/* start in 47:1. Hence 65:2 is 74 */
fygn3 = fygn3[74:73+rows(x)];
fygl = fygl[74:73+rows(x)];
cyr = cyr[74:73+rows(x)];
cmo = cmo[74:73+rows(x)];
/* check that cyr[1] cmo[1] = 65:02 */
sprd = fygn3-fygl;

if levels;

let lbl =  "Y"  "P"  "PCOM"  "NBR"  "TR" "sprd"  "FF"  "M1"  "M2";
/*          1    2     3       4      5     6     7     8     9  */

/* take logs of y,p,nbr,tr,m1,m2 */
x = x[.,1:5]~sprd~x[.,6:8];
x[.,1|2|4|5|8|9] = 100*ln(x[.,1|2|4|5|8|9]);

gosub dovar;

graphjc;
title("Responses to ff shocks");
_plegstr = "Y\000P\000FF";
_psymsiz = 1.5;
_plctrl =   1;
_plegctl = 1;
xy(seqa(1,1,rows(ir)),ir[.,6*N+1|2|7]);

_plegstr = "Y\000P\000PCOM\000NBR\000TR\000sprd\000FF\000M1\000M2";
_plctrl =   1|   1|       0|    0|    0|    0   |   1|    0|    0;
xy(seqa(1,1,rows(ir)),ir[.,6*N+1:6*N+N]);

endif; /* ends if levels */

if levels == 0;     /* specify everything */
                    /* variables: dy, dp, pcom, nbr-tr , tr, dff, dm1, y+p-m2 */

T = rows(x);
x = (100*ln(y[2:T]./y[1:T-1])) ~
    (100*ln(p[2:T]./p[1:T-1]))~
    pcom[2:T]~
    (100*ln(nbr[2:T])-100*ln(tr[2:T]))~
    (100*ln(tr[2:T]))~
    @(ff[2:T]-ff[1:T-1])~@
    sprd[2:T]~
    ff[2:T]~
    (100*ln(m1[2:T])-100*ln(m1[1:T-1]))~
    (100*ln(y[2:T])+100*ln(p[2:T])-100*ln(m2[2:T]));

let lbl =  "DY"  "DP"  "PCOM"  "NBR-TR"  "DTR" "sprd" "FF"  "DM1"  "y+p-m2";

gosub dovar;

/* now unwind cointegration etc to produce ir of variables of interest */

/* variables: dy, dp, pcom, nbr-tr , tr, dff, dm1, y+p-m2 */
ir[.,seqa(1,N,N)] = cumsumc(ir[.,seqa(1,N,N)]);     /* y to log levels */
ir[.,seqa(2,N,N)] = cumsumc(ir[.,seqa(2,N,N)]);     /* p to log levels */
ir[.,seqa(4,N,N)] = ir[.,seqa(5,N,N)]+ir[.,seqa(4,N,N)]; /* nbr = nbr-tr+tr */
@ir[.,seqa(6,N,N)] = cumsumc(ir[.,seqa(6,N,N)]);     /* ff to levels */@
ir[.,seqa(8,N,N)] = cumsumc(ir[.,seqa(8,N,N)]);     /* m1 to log levels */
ir[.,seqa(9,N,N)] = -ir[.,seqa(9,N,N)]+ir[.,seqa(1,N,N)]+ir[.,seqa(2,N,N)];


graphjc;
title("Responses to ff shocks");
_plegstr = "Y\000P\000PCOM\000NBR\000TR\000FF\000M1\000M2";
_plegstr = "Y\000P\000FF";
_psymsiz = 3;
_plctrl =   1|   1|       0|    0|    0|    1|    0|    0;
_plctrl =   1|   1|    1;
_plegctl = 1;
xy(seqa(1,1,rows(ir)),ir[.,5*N+(1|2|7)]);

endif; /* ends if levels */

end;

dovar:
/* subroutine to run vars. Evenutally make a proc */

T = rows(x);
N = cols(x);

lhv = x[k+1:T,.];
rhv = x[k:T-1,.];
                /* rhv is organized x1(t)...xn(t) x1(t-1)....xn(t-1)....*/
i = 2;
do while i <= k;
    rhv = rhv~x[k+1-i:T-i,.];
    i = i+1;
endo;

lhv = lhv -meanc(lhv)';
rhv = rhv - meanc(rhv)';

b = lhv/rhv;

/* b is a matrix of regression coefficients N wide and NK long.
   the left hand variable is on the column, the right hand variable on the row*/

/* give regressions */

format /RD 6,2;
let lbl2 = "lag" ;
j = 1;
do while j <= N;
    "";
    $lbl[j] "regression coefficients";
    $lbl2 $lbl';
    i = 1;
    do while i <= k;
        i b[N*(i-1)+1:N*i,j]';
        i = i+1;
    endo;
    j = j+1;
endo;



err = lhv - rhv*b;
sigma = err'err/T;
q = chol(sigma);

/* q is ut and q'q = sigma */
/*  eps = q'nu(t) E(nunu') = I  expresses the system in orthogonalized form */

/* form vector AR(1) representation */
/* x1(t)                            x1(t-1)           Q'
    ..             ----- b' ----      ...
   xn(t)
   x1(t-1)    =                                       0   nu t
    ..             I    0  0  0                +
   xn(t-1)              I  0  0                       0
    ..
    ..                                                0
   xn(t-k)                 I  0     xn(t-k)


     x(t)    =          A            x(t-1)   +     C  nu(t)
*/

if k > 1;
    A = B'|(eye(N*(k-1))~zeros(N*(K-1),N));
    C = Q'|zeros(N*(K-1),N);
else;
    A = B';
    C = Q';
endif;
@C = eye(N)|zeros(N*(K-1),N);    @


/* Simulate impulse response function */

j = 1;     /* from index */
do while j <= N;
    nu = zeros(N,1);
    nu[j] = 1;
    sim = C*nu;
    simhist = sim[1:N]';
    i = 2;
    do while i <= maxsim;
        sim = A*sim;
        simhist = simhist|sim[1:N]';
        i = i+1;
    endo;
    if j == 1;
        ir = simhist;
    else;
        ir = ir~simhist;
    endif;
    j = j+1;
endo;

/* organizatio of ir
1    2   3    N  N+1 N+2   N+N  2N+1...   2N+N
from 1           from 2         from 3 ....
to1 to2 to3 ...  to1 to2 to3 .. to1 to2 to3 ...
*/
 return;