/* -------------------------------------------------------------------------- */
/*                                 VAR.G
    general purpose proc for doing var's of cointegrated systems in error
    correction form.

    {b,sigma,r,dxir,xir,statcom,e} = var(x,a,lags,horiz,names,bigres,prnt)

    Inputs:
    x     = series, from 0 to T, in levels.  shocks are orthogonalized in the
             order given here.
    a     = cointegrating vector. x*a is
             stationary and is added to the right hand side of the VAR
    lags  = number of lags of differences to include
    horiz = how far out to do impulse response function
    names = vector of names of the cointegrating vectors and then variables
            no more than 3 letters as (t-x) are added in the program.
            e.g. names = Y-C, y, c.
    bigres = 1 extra stuff; 0= bare bones for monte carlo.
    print = 1 to print out results, 0 to not print.

    Outputs:
    b: OLS regression coefficients:
                         dx1(t)   dx2(t)  dx3(t) ...
             const
             ci1(t-1)
             ci2(t-1)
             ..
             dx1(t-1)
             dx1(t-2)
             ..
             dx1(t-lags)
             dx2(t-1)
             ...

     r: choleski decomposition of sigma; sigma = r'r
     dxir: impulse response of dx.
                          response of
                dx1 dx2 dx3 ..    dx1 dx2 dx3 ...
                        to unit shock in
                dx1 dx1 dx1       dx2 dx2 dx2 ...
          t
          t+1
          t+2
          ...
          t+horiz

      statcom: Beveridge Nelson stationary component
      e: error series



    The proc runs a var of the form

    A(L) Dx(t) =  B * (x(t-1)*a) + e.

    by OLS. It uses the choleski decomposition

    var(ee') = sigma = rr'

    so

    e = rv  E(vv') = I

    the final form of the VAR is then

       -1                 -1
    (r)  A(L) Dx(t) =  (r)   B * (x(t-1)*a) + v

    It finds impluse responses by simulating the OLS system (not the constants)
    starting from 0.

*/
/* ------------------------------------------------------------------------- */


proc(7) = var(x,a,lags,horiz,names,bigres,prnt);
local dx,t,rhv,i,j,b,e,sigma,r,seb,R2,AR,offst,G,mu,rinv, Arorth,Gorth,muorth,
       nseries,blank,header,table,mask,fmt,d,adder,xl,
       dxsim, xsim, compan, vshock, shock, dxir, xir,
       Amat, arows, ll, lr, rl, rr, bnocon, cofst, statcom,
       fstat, pvf, lastrh,
       ci, dci, lhv, bci, eci, sigci, sebci, tci, ss, correl, vardec, denom,
        dfcon,rhotest;

     statcom = 0;

/* ------ */
/* checks */
/* ------ */

    if rows(a) /= cols(x);
        "proc var: cointegrating vector not conformable to data";
    endif;
    if rows(names) /= cols(a)+cols(x);
        "proc var: names not conformable to data, ci vectors";
    endif;

/* ------------------------------------------------------------ */
/*             set up right hand variables                      */
/* rhv = const ci(t-1) .. x1(t-1) x1(t-2) .. x2(t-1) x2(t-2) .. */
/* ------------------------------------------------------------ */

    Dx = x[2:rows(x),.]-x[1:rows(x)-1,.];
    xl  = x[1:rows(x)-1,.];                    /* note x is LAGGED levels */
    T = rows(xl);
    nseries = cols(xl);


    Rhv = ones(T-lags,1)~(xl[1+lags:T,.]*a);    /* constant, lagged ci vec*/
    if lags > 0;                                /* then lags of dx */
        i = 1;
        do while i <= cols(xl);
            j = 1;
            do while j <= lags;
                Rhv = Rhv~dx[1+lags-j:T-j,i];
                j = j+1;
            endo;
            i = i+1;
         endo;
     endif;

/* ----------------------------- */
/* run regression of dx's on rhv */
/* ----------------------------- */

   b = dx[1+lags:T,.]/rhv;                 /* b is a matrix.  each column gives
                                               ols of that lhv on all rhv's  */
   e = dx[1+lags:T,.]-rhv*b;
   sigma = e'e/(T-lags);
   seb = (diag(inv(rhv'rhv))*diag(sigma)')^(1/2);
   R2 = 1 - diag(sigma)./(stdc(dx[1+lags:T,.])^2);
   fstat = (r2./(1-r2)).*(T-cols(xl)-1)/cols(xl);
   pvf = 0*fstat;
   i = 1;
   do while i <= rows(fstat);
        pvf[i] = cdffc(fstat[i],cols(xl),T-cols(xl)-1);
        i = i+1;
   endo;


/* --------------------------------- */
/*  A(L)Dx(t) = mu + G(a'x(t-1)) + e */
/*                representation.    */
/*  AR is A0~A1~A2...                */
/* --------------------------------- */
  if bigres;
      AR = eye(cols(dx));
      offst = cols(a) + 1;
      i = 1;
      do while i <= lags;
        j = 1;
        do while j <= cols(dx);
            AR = AR ~ (-b[offst+i+(j-1)*lags,.]');
            j = j+1;
        endo;
        i = i+1;
      endo;

      G =  b[2:cols(a)+1,.]';
      mu = b[1,.]';
  endif;
/* ----------------------------- */
/*  orthogonalize in order given */
/* ----------------------------- */

   r = chol(sigma)';
   rinv = inv(r);

   if bigres;
       arorth = rinv*ar;
       gorth = rinv*g;
       muorth = rinv*mu;
    endif;

/* -------------------------------------- */
/* find impulse response function         */
/* this uses the nonorthogonalized        */
/* representation, finds the e innovation */
/* corresponding to a unit v innov.       */
/* -------------------------------------- */

   /* ---------------------------------------- */
   /* simulate ols system in companion form   */
   /* starting from sample means.             */
   /* simulation ignores means.               */
   /* [ Dy(t) ] =  [ b' ] [ y(t-1) - c(t-1) ]   +  [ e1 ]
      [ Dc(t) ]           [ Dx(t-1)         ]      [ e2 ]
                          [ Dx(t-lag)       ]
                          [ Dy(t-1)         ]
                          [ Dy(t-lag)       ]

      xsim              compan

      start from y = c = 0 and all diffs = 0   */
   /* ---------------------------------------- */

   shock = 1;                                    /* which series gets shock? */
   do while shock <= nseries;
        gosub dosim;
        if shock == 1;
            dxir = dxsim;
            xir = xsim;
        else;
            dxir = dxir~dxsim;
            xir = xir~xsim;
         endif;
        shock = shock+1;
   endo;

   vardec = recserar(dxir^2,dxir[1,.]^2,ones(1,4));
                         /* this is decomposition of forecast error variance */

/* ------------------------------------ */
/* find beveridge nelson decomposition  */
/* again, you can simply ignore means   */
/* ------------------------------------ */

/* form companion matrix: updates rhv(t)' = A*rhv(t-1)'+ error */
  if bigres;
    statcom = 0;
    if cols(x) == 2;  /* not yet implemented for more right hand vars */

        arows = rows(b)-1;
        bnocon = b[2:rows(b),.];
        Amat = zeros(arows,arows);
        Amat[1,.] =  a'(bnocon') + (1~zeros(1,arows-1));     /* ONE  ci vect */
        if lags >= 2;
            j = 1;
            do while j <= nseries;
            Amat[(2+(j-1)*lags),.] =  bnocon[.,j]';           /* Dxj(t)       */
                ll = (3+(j-1)*lags) ;
                lr = (3+(j-1)*lags+lags-2);
                rl = (2+(j-1)*lags);
                rr = (2+(j-1)*lags+lags-2);
                Amat[ll:lr,rl:rr] = eye(lags-1);
                j = j+1;
            endo;
            cofst = zeros(cols(a)+lags,1)|1|zeros(lags-1,1);

        endif;

        if lags == 1;
            Amat[2:arows,.] = bnocon';
            cofst = zeros(cols(a)+1,1)|1;
        endif;

        cofst = cofst'Amat*inv(eye(arows)-amat);
        cofst = 0~cofst;
                                                     /* rhv(t) is t-1 info ! */
        statcom = -(rhv[2:rows(rhv),.]-meanc(rhv)')*cofst';
        lastrh = 1~(x[rows(x),.]*a);
        if lags > 0;                                /* then lags of dx */
            i = 1;
            do while i <= cols(xl);
                j = 1;
                do while j <= lags;
                    lastRh = lastRh~dx[T+1-j,i];
                    j = j+1;
                endo;
                i = i+1;
            endo;
        endif;

        statcom = statcom |
                  (-(lastrh - meanc(rhv)')*cofst');

    endif;
  endif;

/* -------------------------------- */
/* augmented Dickey-Fuller tests on */
/* ci vector                        */
/* -------------------------------- */

    if bigres;
        "";
        "Augmented Dickey fuller regressions of c-y";
        "D(c-y)(t) = const. + rho(c-y)(t-1) + b1 D (C-y) (t-1) + ... ";
        format /RD 9,3;

        ci = xl*a;
        dci = ci[2:T]-ci[1:T-1];      /* ci is lagged rel. to dci */

        j = 0;
        do while j <= 5; /* max amount of augmenting */
             " ";
            lhv = dci[1+j:T-1];
            rhv = ones(T-1-j,1)~ci[1+j:T-1];
            if j > 0;
                i = 1;
                do while i <= j;
                    rhv = rhv~dci[1+j-i:T-i-1];
                i = i+1;
                endo;
            endif;
            bci = lhv/rhv;
            eci = lhv-rhv*bci;
            sigci = eci'eci/(T-j-1);
            sebci = (diag(inv(rhv'rhv))*diag(sigci)')^(1/2);
            tci = bci./sebci;

            if j == 0;
                dfcon = 1;
            endif;
            if j > 0;
                dfcon = 1-sumc(bci[3:rows(bci)]);
            endif;
            rhotest = bci[2]/dfcon*(T-1);

            "coeff. " bci';
            "se     " sebci';
            "tstat  " tci';
            "rhostat" rhotest;
            j = j+1;
        endo;
    endif;

/* ------------------------------ */
/*    print table of results      */
/*  some have been commented out  */
/* ------------------------------ */

 if prnt;
   "";
   "VAR estimates";
   "";

    /* ----------------------------------- */
    /* create header to identify variables */
    /* ----------------------------------- */

    let blank = "(t-1)" "(t-2)" "(t-3)" "(t-4)" "(t-5)" "(t-6)" "(t-7)" "(t-8)";
    let header = "    " "const" ;
    j = 1;
    do while j <= cols(a);
        header = header|(names[j] $+ "(t-1)");     /* add names of ci vectors */
        j = j+1;
    endo;

    i = cols(a) + 1;
    do while i <= rows(names);
        j = 1;
        do while j <= lags;
            header = header|(names[i] $+ blank[j]);
            j = j+1;
        endo;
        i = i+1;
    endo;

    /* -------------------------- */
    /* create table of ols coeffs */
    /* -------------------------- */

    table = header'|
            ((names[cols(a)+1:rows(names)] $+ "(t)") ~ b');
    fmt = "lf"~9~3;
    mask = zeros(1,cols(table))|(zeros(rows(b'),1)~ones(rows(b'),cols(b')));

    "OLS coefficients";
    call printfm(table, mask, fmt);
    " ";

    "Standard Errors";
    table = header'|
            ((names[cols(a)+1:rows(names)] $+ "(t)") ~ seb');
    call printfm(table, mask, fmt);
    "";
    "t Ratios";
    table = header'|
            ((names[cols(a)+1:rows(names)] $+ "(t)") ~ (b'./seb'));
    call printfm(table, mask, fmt);
    "";

    "R2, pvalue";
    format /RD 10,4; R2~pvf;

    "Variance/covariance, correlation matrix of errors";

    ss = diag(sigma)^(1/2);
    correl = sigma ./ (ss*ss');
    "sigma" sigma;
    "corr " correl;

    /* ------------------------------------ */
    /* AR and orthogonalized representation */
    /* ------------------------------------ */
    if bigres;
        "A(L)Dx(t) = mu + G(a'x(t-1)) + e representation. ";
        "A(L): A0 ~ A1 ~ A2 ...";
        format /RD 9,3;
        AR;
        "G";G;
        "mu";mu;

        "A(L)Dx(t) = mu + G(a'x(t-1)) + v orthogonalized representation. ";
        "This is the same as ols with current consumption in y ";
        "A(L): A0 ~ A1 ~ A2 ...";
        format /RD 9,3;
        ARorth;
        "G";Gorth;
        "mu";muorth;
    endif;

    /* ------------------------------------------------ */
    /* print variance-covariance and correlation matrix */
    /* ------------------------------------------------ */

    "";
    "Variance covariance matrix";
    format /RD 9,2;
    sigma;
    "";
    "Correlation matrix";
    d = diag(sigma);
    sigma./((d*d')^(1/2));
    "";
    "choleski decomp of sigma";
    r;
    "";

    /* ---------------------------------- */
    /* print impulse-response functions   */
    /* ---------------------------------- */
/*
    "Impulse response function";
    "Differences";
    "   x1 to x1 shock, x2 to x1 shock.. x1 to x2 shock x2 to x2 shock..";
    "Differences: Dx1 Dx2.. to shock in x1,x2,...";
    format /RD 10,4; dxir;
    "";
    "Levels";
    "   x1 to x1 shock, x2 to x1 shock.. x1 to x2 shock x2 to x2 shock.. ";
    format /RD 10,4; xir;
*/

    "Decomposition of variance";
    "order:x1 to x1 shock, x2 to x1 shock ... x1 to x2 shock x2 to x2 shock.." ;
    "1: one step ahead forecast error variance";
    vardec[1,.];
    "variance ratio to series variance x 100";
    denom = (1~1) .*. stdc(dx)';
    100*vardec[1,.]./(denom^2);
    "";
    "2: horiz step ahead forecast error = series variance";
    vardec[rows(vardec),.];
    "variance ratio to series variance x 100";
    100*vardec[rows(vardec),.]./(denom^2);
    "";
    "series variance";
    denom^2;
    "series standard deviation";
    denom;


  endif;   /* ends if on print section */

retp(b,sigma,r,dxir,xir,statcom,e);

/* ------------------------------------- */
/* dosim: calculates responses to shocks */
/* ------------------------------------- */

dosim:
   dxsim = zeros(horiz,nseries);                   /* setup: differences     */
   xsim = dxsim;                                   /* levels                 */

   /* t=1 levels and diffs  */

   compan = zeros(cols(a)+lags*nseries,1);         /* last period's state    */
   vshock = zeros(nseries,1);
   vshock[shock] = 1;                              /* unit shock             */
   dxsim[1,.] =  (r*vshock)';                      /* e shock equal to v shk */
   xsim[1,.] =  dxsim[1,.];

   /* ---------------------- */
   /* simulate resp to shock */
   /* ---------------------- */

   i = 2;
   do while i <= horiz;

      /* ------------------------------------------ */
      /* update compan (compan is conformable to b) */
      /* ------------------------------------------ */

      compan[2+cols(a):rows(compan)] =               /* shift down */
                 compan[1+cols(a):rows(compan)-1];
      compan[1:cols(a)] = (a'(xsim[i-1,.]')) ;       /* new lagged ci vector  */
      j = 1;                                         /* new lagged dx(t)'s    */
      do while j <= nseries;
        compan[1+cols(a)+(j-1)*lags] = dxsim[i-1,j];
        j = j+1;
      endo;

      /* ------------- */
      /* generate Dx,x */
      /* ------------- */

      dxsim[i,.] = (b[2:rows(b),.]'compan)' ;
      xsim[i,.]  = xsim[i-1,.] +  dxsim[i,.];

      i = i+1;
   endo;
return;

endp;


/* test call

output file = var.out reset;

load x[171,11] = gene.dat;
 GCN82   =x[.,4];
 GCS82   =x[.,5];
 GNP82   =x[.,10];
 x = 0;

 gcns82 = gcn82+gcs82;

 lags = 2; horiz = 50;
 x = 100*ln(gcns82~gnp82);
 a = 1|-1;
 names = "C-Y"|"c"|"y";
 prnt = 1;
 {b,r,dxir,xir,statcom} = var(x,a,lags,horiz,names,prnt);

output off;   */