/* -------------------------------------------------------------------------- */
/*                                 VAR2.G
    copied from var.g.  does things with current consumption in GNP equation.

    {b,sigma} = var2(x,a,lags,horiz,names,prnt)

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


proc(2) = var2(x,a,lags,horiz,names,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,
       rhv1,
       bc, ec, sigmac, sebc, R2c, fstatc, pvfc,
       by, ey, sigmay, seby, R2y, fstaty, pvfy;


/* ------ */
/* 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) .. */
/* rhv1 = const ci(t-1) .. x1(t) 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);


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

     rhv = rhv1[.,1:2]~rhv1[.,4:cols(rhv1)];   /* with no contempo cons */

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

   bc = dx[1+lags:T,1]/rhv;
   ec = dx[1+lags:T,1]-rhv*bc;
   sigmac = ec'ec/(T-lags);
   sebc = (diag(inv(rhv'rhv))*diag(sigmac)')^(1/2);
   R2c = 1 - diag(sigmac)./(stdc(dx[1+lags:T,1])^2);
   fstatc = (r2c./(1-r2c)).*(T-cols(xl)-1)/cols(xl);
   pvfc = cdffc(fstatc,cols(xl),T-cols(xl)-1);

   by = dx[1+lags:T,2]/rhv1;
   ey = dx[1+lags:T,2]-rhv1*by;
   sigmay = ey'ey/(T-lags);
   seby = (diag(inv(rhv1'rhv1))*diag(sigmay)')^(1/2);
   R2y = 1 - diag(sigmay)./(stdc(dx[1+lags:T,2])^2);
   fstaty = (r2y./(1-r2y)).*(T-cols(xl)-1)/cols(xl);
   pvfy = cdffc(fstaty,cols(xl),T-cols(xl)-1);

   b = (bc[1:2]|0|bc[3:rows(bc)])~by;
   seb = (sebc[1:2]|999|sebc[3:rows(bc)])~seby;
   e = ec~ey;
   sigma = (e'e)/(T-lags);

/* ------------------------------------- */
/* print stuff moved to here temporarily */
/* ------------------------------------- */
if prnt;
   "";
   "VAR estimates with current consumption in GNP equation";
   "These are only for a potential table, use VAR.G for all calculations";
   "Note that this produces the same numbers as teh orthogonalized var";
   "with OLS standard errors";
   "";

    /* ----------------------------------- */
    /* 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";
    header = header | (names[1] $+ "(t-1)") | (names[2] $+ "(t)");

    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);
    "";

    format /RD 10,4;
    "R2, c, y" R2c r2y;
    "pvalue c,y" pvfc pvfy;

    "Variance/covariance, correlation matrix of errors";

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

endif;   /* ends print if */

retp(b,sigma);

endp;