/*******************   VAR code    ***************************
Author: Paul Fackler <FACKLER@econbus2.econ.ncsu.edu>
Date:  7 Mar 1996
Provided free, without guarantees, for public non-commercial use.

There have been several postings lately about VAR methods.  It seems
that there are several out there but some of these are not provided
as source code and some are not free.  Here some fairly
plain code to compute VARS, IRFs, and FEVDs that I wrote a number of
years ago.  It gets the job done and is open to examination.  It is
modular and allows for alternative identification schemes besides the
usual Cholesky factorization.  I have other code to compute maximum
likelihood estimates for non-triangular identifications and will be
glad to send or post it if there is interest. (I published two papers
on this topic in the 80s and would be happy to rid myself of some
old reprints on a first come first served basis.)

Description:
  Procedure for doing vector autoregressions (VARs).  This procedure
  calculates an AR representation of a vector time series, i.e. a
  regression of a set of variables on a specified number of lags of
  all the variables.

  INPUTS:
    X : an NxK data matrix. This includes just the endogenous
variables
        of the model.
    LAGS : the number of lags desired in the AR representation.  The
first
        LAGS observations are used only for the regressors, hence
there
        will effectively be N-LAGS observations in this problem.
    Z : an NxP or (N-LAGS)xP matrix of deterministic (or exogenous)
        variables (these might include a constant, time trends, or
        dummy variables as well as stochastic but exogenous
variables).
        If there are N rows in the matrix, the first LAGS rows will be
        ignored.  Other legal values for Z are 0 and 1.  If Z=0 then
        no deternimistic variables will be used.  If Z=1 then a
constant
        term will be included.
    NAMEOUT : a 6 letter or less prefix used for names of GAUSS matrix
        files saved by the program.

  OUTPUT
    AR: the regression coefficients, stored as a (K*LAGS+P)xK matrix.
        This may be thought of as LAGS KxK matrices, the ith of which
is
        associated with the ith lag, with a PxK matrix concatenated
at its
        bottom.
        Note that these matrices postmultiply the variables so
columns are
        associated with equations, rows with regressors.

  SAVED TO DISK:
    XX: the inverse moment matrix (inv(moment([X(t-1)...X(t-LAGS)
Z]))).
        This is useful in calculating the sampling variances of the
        coefficients when multipied by the diagonal elements of CV.
        If the problem is large (# of regressors>90) then only the
        diagonal elements of XX are saved.
    CV: the sample covariance matrix of the errors, a KxK matrix.
Note that
        this is the ML estimator, calculated as U'U/(N-LAGS).
    U : the regression residuals, an (N-LAGS)xK matrix.
  Each of these matrices will be saved to disk with the prefix sent
  followed by the names given above.  Thus if the prefix sent is "R1"
  then the following matrices will be stored on disk: R1XX, R1CV,
  R1U. These may be accessed by using the LOAD command.


*************************************************************************/
proc var(x,lags,z,nameout);
  local k,kl,i,j,ik,jk,n,numdet,numreg,xx,xy,ar,u,vcv,name,
    l1,l2,kl1,kl2,xx1,xx2,xx3,xy1,xy2,temp,d1,d2,dyy,dzz;
  cls;
@
 Calculates problem dimensions and checks for errors in problem size,
 data, or in room to write output disk files
@
  k=cols(x);
  n=rows(x);
  if rows(packr(x))/=n;
    "Error: Missing values encountered in X data";
  stop;endif;
  if rows(z)>1;
    if rows(z)/=n-lags;
      if rows(z)==n;z=z[lags+1:n,.];
      else;
      "Error: Rows(Z) must equal rows(X)-lags";
    stop;endif;endif;
    if rows(z)/=rows(packr(z));
      "Error: Missing values encountered in deterministic data (Z)";
    stop;endif;
    numdet=cols(z);
  elseif z==0;
    numdet=0;
  else;
    z=ones(n-lags,1);
    numdet=1;
  endif;
  lags=trunc(lags[1,1]);
  if lags<1;"Error: Improper lag specification";stop;endif;
  kl=k*lags;
  numreg=kl+numdet;
@  if numreg>180;"Error: Problem too big";stop;endif;@
  nameout=upper(nameout);
  i=vals(nameout);
  if rows(i)>2;
    if i[2,1]==58;i=i[1,1]-64;else;i=0;endif;
    if i;i=2;endif;
    if strlen(nameout)>6+i;
      nameout=strsect(nameout,1,6+i);
      "Note: Output file name prefix shortened to " $ nameout;
    endif;
  endif;
  dyy=sqrt(diag(moment(trimr(x,lags,0),0))');
  x=x./dyy;
  if numdet;
    dzz=sqrt(diag(moment(z,0))');
    z=z./dzz;
  endif;

@
 Calculates the coefficient matrix, AR, a (numreg x k) matrix.
 This matrix is written to disk.
 Inv(x'x), a (numreg x numreg) matrix, stored as xx, is written to
disk.
 XX and XY are cleared to open work space for further calculations.
 Covariance matrix of the ith column of AR calculated by multipling
 xx by the variance of the residuals of the ith equation.
@
    xx=zeros(numreg,numreg);xy=zeros(numreg,k);
    i=0;do until i==lags;i=i+1;
      ik=i*k;
      xy[ik-k+1:ik,.]=trimr(x,lags-i,i)'trimr(x,lags,0);
      if numdet;
        xx[ik-k+1:ik,kl+1:numreg]=trimr(x,lags-i,i)'z;
        xx[kl+1:numreg,i*k-k+1:i*k]=xx[ik-k+1:ik,kl+1:numreg]';
      endif;
      j=i-1;do until j==lags;j=j+1;
        jk=j*k;
        xx[ik-k+1:ik,jk-k+1:jk]=trimr(x,lags-i,i)'trimr(x,lags-j,j);
        if j/=i;
          xx[jk-k+1:jk,ik-k+1:ik]=xx[ik-k+1:ik,jk-k+1:jk]';
        endif;
      endo;
    endo;
    if numdet;
      xx[kl+1:numreg,kl+1:numreg]=moment(z,0);
      xy[kl+1:numreg,.]=z'trimr(x,lags,0);
    endif;
    name=nameout$+"XX";save ^name=xx;
    name=nameout$+"XY";save ^name=xy;
    xx=invpd(xx);
    ar=xx*xy;xy=0;
 @
  Calculates the residuals, stored in u.
    These are written to nameu and cleared.
  Residual covariance matrix stored as VCV. These are written to
namecv.
 @
  u=trimr(x,lags,0);
  if numdet;
    u=u-z*ar[kl+1:numreg,.];
  endif;
  i=0;do until i==lags;i=i+1;
    u=u-trimr(x,lags-i,i)*ar[(i-1)*k+1:i*k,.];
  endo;
  clear x,z;
  u=u.*dyy;
  vcv=moment(u,0)/(n-lags);
  name=nameout$+"U";save ^name=u;
  name=nameout$+"CV";save ^name=vcv;
  clear u,vcv;
  ar=ar.*dyy;
  ar[1:kl,.]=ar[1:kl,.]./(ones(lags,1).*.dyy');
  if numdet; ar[kl+1:numreg,.]=ar[kl+1:numreg,.]./dzz'; endif;
  retp(ar);
endp;

@
  A procedure for calculating a vector moving average representation
from
  a vector autoregressive representation.

INPUTS
    AR: the regression coefficients, stored as a (K*LAGS+P)xK matrix.
       These may be calculated by PROC VAR.
    LAGS : the number of lags desired in the AR representation.  The
first
        LAGS observations are used only for the regressors, hence
there
        will effectively be N-LAGS observations in this problem.
    MALAGS : the number of lags desired in the MA representation and
        hence in the impulse response function.

OUTPUT
    MA: the MA representation associated with the innovations, i.e.
the
        step ahead forecast errors.
  Note this is stored as an (MALAGS*KxK) matrix. The ith block of k
rows
  is the matrix of coefficients associated with the ith lagged
innovation.
  Also note the convention followed here is that u (the innovation
vector)
  is a row vector and the coefficients therefore postmultipy u.

@
proc vma(ar,lags,malags);
  local i,j,k,ik,ma;
  k=cols(ar);
  ma=zeros(malags*k,k);
  i=k*minc(lags|malags);
  ma[1:i,.]=ar[1:i,.];
  i=0;do until i==malags-1;i=i+1;
    j=k*minc((i+lags)|malags); ik=i*k;
    ma[ik+1:j,.]=ma[ik+1:j,.] + ar[1:(j-ik),.]*ma[ik-k+1:ik,.];
  endo;
  retp(ma);
endp;

@
  Computes the A0 matrix for a recursive identification.
  This may be passed to IMPULSE to obtain the IRF associated with
  such a model, etc.  The Cholesky decompostion of the innovation VCV
  matrix is used.
  INPUTS:
    VCV: the innovation covariance matrix (called SIGMA in VAR2)
    IND: a K-vector listing the ordering of the recursive model
  OUTPUTS:
    A0: the contemporaneous coefficients matrix (inv(A0*A0')=VCV)
@
proc reca0(vcv,ind);
  local a0;
  vcv=vcv[ind,ind];
  a0=vcv;
  a0[ind,ind]=inv(chol(vcv));
  retp(a0);
endp;
@
  Computes the Impulse Response Function (IRF) of a K variable VAR
  with an orthogonalization defined by A0, the matrix of contemporaneous
  coefficients. Input is the Moving Average (MA) representation of the
  system and the A0 matrix.  The former can be obtained from the VAR
  procedure, while the latter is a transformation of the VAR residual
  covariance matrix (VCV), which is placed in memory by the VAR
procedure.
  A general set of routines for calculating A is available in the file
  VAR2, which also uses PROC MINMUM and associated procedures. (Note
VCV
  is called SIGMA in those routines).
  RECA0 can be used to obtain a triangular A matrix.
  This proc returns the IRF as an (MALAGS+1)x(K*K) matrix, made up of
  K (L+1)xK blocks, the ith of which gives the response of the the K
  variables to a unit shock in the ith impulse.  Note that
  row i of the IRF correspond to the responses to shocks at a lag of
  i-1.  Hence the 1st row is the contemporaneous response.
  A standard output is printed if OUTCODE is non-zero.
@
proc impulse(ma,a0,outcode);
  local i,j,k,l,m,s;
  cls;
  k=rows(a0);
  if cols(a0)/=k;"Error: A0 is not square";stop;endif;
  l=rows(ma)/k;                           @ This is MALAGS of PROC
VAR @
  if l/=trunc(l) or cols(ma)/=k;
    "Error: MA and A0 are not conformable";stop;endif;
  a0=inv(a0);
  i=0;do until i==l;i=i+1;
    ma[(i-1)*k+1:i*k,.]=a0*ma[(i-1)*k+1:i*k,.];
  endo;
  ma=reshape((a0|ma),l+1,k*k);

  if outcode;
  s=seqa(1,1,k);
  format 8,0;
  "\f";
  "                         IMPULSE RESPONSE FUNCTIONS";
  "                            response of variable:";
  "  to    at";
  "shock  lag" s';
  "-------------";;
  j=0;do until j==k;j=j+1;"---------";;endo;"";
  j=0; do until j==k;j=j+1;
    i=0;do until i==l+1;
      format 4,0;
      if i==0;j;;" " i;;
      else;"      " i;;endif;
      format 8,4;"  ";;ma[i+1,s];
      i=i+1;
    endo;
    "\f";
    s=s+k;
  endo;
  endif;
  retp(ma);
endp;

@
  Computes the Forecast Error Decomposition (FED) for orthogonal
shocks.
  Input is the Impulse Response function for the system.  This can
  be obtained from the VAR procedure or from the IMPULSE procedure.
  This proc returns the FED as an (MALAGS+1)x(K*K) matrix.  This is K
  (MALAGS+1)xK blocks, the ith of which gives the decomposition for
the
  ith variable.  The ith row of each block is associated with the i-
step
  ahead forecast, the jth column is associated with the jth impulse.
  A standard output is printed if OUTCODE is non-zero.
@
proc fedecomp(irf,outcode);
  local i,j,d,k,kk,l,t,q,v,s,var;
  cls;
  l=rows(irf);                              @ Equals MALAGS+1 @
  k=round(sqrt(cols(irf)));
  kk=k*k;
  irf=reshape(irf.*irf,k*l,k);
  d=zeros(k,k);v=zeros(l,kk);
  var=zeros(l,k);
  t=0;do until t==l;t=t+1;
    d=d+irf[(t-1)*k+1:t*k,.];
    var[t,.]=sumc(d)';
    v[t,.]=reshape((100*d./var[t,.])',1,k*k);
  endo;

  if outcode;
  s=seqa(1,1,k);
  format 5,0;
  "                               FORECAST ERROR DECOMPOSITION";
  "                         % of forecast error variance due to
shock:";
  "  for    steps     std.";
  "variable ahead     err.  "  s';
  "------------------------";;
  j=0;do until j==k+1;j=j+1;"------";;endo;"";
  j=0; do until j==k;j=j+1;
    i=0;do until i==l;i=i+1;
      format 4,0;
      if i==1;j;;"   " i;;
      else;"        " i;;endif;
      format 8,4;"    ";;sqrt(var[i,j]);;
      format 5,1; v[i,s];
    endo;
    "";
    s=s+k;
  endo;
  endif;

  retp(v);
endp;

@ This procedure calculates 1 through m step ahead projections
starting at
  period t (i,e, the tth observation of X, the data matrix).  AR are
the
  autoregressive coefficients calculated by the VAR procedure.  Z is a
  MxP matrix of deterministic variables corresponding to the
projection
  period.  If Z=0 then no deterministic variables are used. If Z=1
then
  an m vector of ones is used. If Z and X have the same number of
rows then
  rows (t+1) through (t+m) of Z will be used. Other than these cases
  Z must be MxP or an error message will be printed.
  Output is an MxK matrix with the i,jth element representing the
  expectation at time t of the jth variable at time t+i.
@
proc proj(x,ar,t,m,z);
  local i,ey,k,l,kl,temp,r;
  if z==1;z=ones(m,1);endif;
  if rows(z)==rows(x);z=z[t+1:t+m,.];endif;
  if t>rows(x);"Error: Starting observation,t, greater than rows(X)";
     stop; endif;
  if rows(z)/=m and z/=0;
     "Error: Deterministic variables should contain m
observations";stop;endif;
  k=cols(ar);
  if z==0;l=rows(ar)/k;else;l=(rows(ar)-cols(z))/k;endif;
  kl=k*l;
  if l/=trunc(l);"Error: cols(Z) doesn't match size of AR";stop;endif;
  ey=rev(x[t-l+1:t,.]);
  i=0;do until i==m;i=i+1;
     if z==0;ey=(reshape(ey[1:l,.],1,kl)*ar)|ey;
     else;ey=((reshape(ey[1:l,.],1,kl)~z[i,.])*ar)|ey;
     endif;
  endo;
  retp(rev(ey[1:m,.]));
endp;

@
  Calculates the roots of an AR model and calculates how many, if
any, lie
  outside the unit circle (indicating non-stationarity).
  INPUTS:
    AR: the AR coefficients matrix obtained from proc VAR
    LAGS: the number of lags in the VAR
  OUTPUTS:
    R: a K*LAGS x 3 matrix containing the real and imaginary parts of
the
       roots of the model and their modulus
  Note: the external function EIG must be previously loaded.
@
proc roottest(ar,lags);
  local r,rv,iv,k,kl;
  k=cols(ar);
  kl=k*lags;
  ar=ar[1:kl,.]~(eye(kl-k)|zeros(k,kl-k));
  {rv,iv}=eigrg(ar);
  r=sqrt(rv.*rv+iv.*iv);
  k=sumc(r.>=1);
  k "root(s) outside the unit circle";
  retp(rv~iv~r);
endp;
