/* Structural VARs:
Essentially this code sets up the likelihood
function, score and information matrix for a VAR that is identified
using contemporaneous coefficient restrictions, which may be thought
of as restrictions on the error covariance matrix.  It is required
that the reduced form (VAR) error covariance matrix is estimated
(e.g., using the VAR file I posted to the American University
Archive). Also matrices defining the appropriate coefficient
restrictions must be specified (only linear restrictions are allowed;
most applications have used zero restrictions).  There are utilities
designed to facilitate the creation of the appropriate restriction
matrices.

I wrote this code many years ago and at that time I arranged matrices
so columns refered to equations and rows to variables (i.e. the
problem was written in the form y(t)A(0)=sum y(t-i)A(i)+e, where y is
a row vector).  I have since transposed my thinking to conform to
standard practice but I don't recall ever rewritting the code.  In
any case, beware when you use this code.

A good reference on this approach is Carlo Giannini, "Topics in
Structural VAR Econometrics" Springer-Verlag Lecture Notes in
Economics and Mathematical Systems #381, 1992.

Good Luck.
*/

/************
SVAR - Structural Vector Autoregression
Coded by Paul L. Fackler, North Carolina State University

This code is provided "as is" with no claims as to its suitablility
for specific estimation tasks.  It may be freely copied and used so
long as its source is appropriately referenced.
*************/

@
  This file defines the necessary procedures to calculate the log-
  likelihood and gradient vector for the optimization routine to
  use. It assumes that  the matrices SIGMA, K, CFACTOR, Z1, Z2, W1,
  and W2 are defined globally (the Z and W makrices can be created
  using ZMAKE).
  SIGMA=U'U/T where u are the residuals from the VAR estimation
  (the innovations) and T is the number of sample obseravations.
  K=rows(SIGMA) is the dimension of the model, while CFACTOR is the
  constant term in the loglikelihood function.  Note that the
  loglikelihood, score and information matrix values computed are all
  divided by T.  This is especially important when computing model
  selection criteria or estimators covariances.  In both case the
  values returned should be multiplied by T.  The matrices Z1, Z2,
  W1, and W2 define the identifying (and over-identifying)
  restrictions on A and B according to vec(A)=Z1*THETA + W1 and
  vec(B)=Z2*THETA + W2, where THETA is a vector of free parameters
  to be estimated.
  These procedures can then be passed to a maximization routine to
  compute the ML estimates.
@
if scalmiss(typecv("sigma"));
  clear sigma,k,z1,z2,w1,w2,cfactor;
endif;
proc lkinfo;retp("");endp;   @ dummy proc definition @
clear a,b;

@
  This proc creates the Z1, Z2, W1, and W2 restriction matrices for
  an identification involving only zero restrictions and
  normalization. IND1 and IND2 specify the elements of VEC(A) and
  VEC(B), respectively, that are NON-ZERO.  If all elements are to
  be zero then pass a zero.
  NORM is one if normalization (ones on the diagonal) is applied to
  A, otherwise B is normalized.  Note that typically none of the
  diagonal elements of the normalized matrix should be included in
  the associated IND.  The proc also checks the model identification
  by passing a random vector to LKINFO.  This returns the information
  matrix, which is checked for posite definiteness.
@
proc zmake(ind1,ind2,norm);
  local p1,p2,kk,n,temp;
  if rows(sigma)/=cols(sigma);
    "Error: SIGMA is not square"; stop;
  endif;
  if k/=rows(sigma);
    "Error: K and SIGMA are not conformable"; stop;
  endif;
  cfactor=0.5*k*ln(2*pi);
  kk=maxc(maxc(ind1)¦maxc(ind2));
  if (k*k)<kk;
    "IND1 or IND2 incompatible with size of K"; stop;
  endif;
  kk=k*k;
  p1=rows(ind1); p2=rows(ind2);
  if ind1/=0;
    z1=eye(kk);
    z1=z1[.,ind1];
    if ind2/=0; z1=z1~zeros(kk,p2); endif;
    n=p1;
  else;
    z1=zeros(1,1); a=eye(k);
    n=0;
  endif;
  if ind2/=0;
   z2=eye(kk);
   z2=z2[.,ind2];
   if ind1/=0;z2=zeros(kk,p1)~z2;endif;
   n=n+p2;
  else;
   z2=zeros(1,1); b=eye(k);
  endif;
  if n==0;
    "Error: Dimension of THETA equals zero"; stop;
  endif;
  if norm; w1=vec(eye(k)); w2=0; else; w2=vec(eye(k)); w1=0; endif;
  trap 1;
  temp=lkinfo(rndn(n,1));
  if scalmiss(temp) or scalerr(invpd(temp))==20;
    "Warning: Check model identification.";
  endif;
  trap 0;
  retp("");
endp;

@
  Returns the negative of the maximum likelihood value
  divided by T, the number of observations.  This is the value that
  the LK procedure is attempting to attain for a given
  (exact) identification.
@
proc lksig;
  k=rows(sigma);
  cfactor=0.5*k*ln(2*pi);
  retp(cfactor + 0.5*(ln(det(sigma))+rows(sigma)));
endp;


@
  Basic procedures for computing the loglikelihood, score and
 information matrix for VAR problems.  It is assumed the SIGMA, K,
  CFACTOR, Z1, Z2, W1, and W2 are correctly defined.  Note that
  these functions should be multiplied by T, the number of
  observations, to be correct.  Alternately they can be thought of
  as averages.
  Note that in each of these three functions there are three
  algorithms for the cases in which B=I, A=I or neither A nor B=I.
  The speed of calculation runs in this order.
@

proc lk(theta);
  local binv,abinv;
  if z2==0;
    a=reshape(z1*theta+w1,k,k)';
    retp(cfactor + 0.5*vec(a)'vec(sigma*a) - ln(abs(det(a))));
  elseif z1==0;
    b=reshape(z2*theta+w2,k,k)';
    binv=inv(b);
    retp(cfactor+ 0.5*vec(binv)'vec(sigma*binv)-ln(abs(det(binv))));
  else;
    a=reshape(z1*theta+w1,k,k)';
    b=reshape(z2*theta+w2,k,k)';
    abinv=a*inv(b);
  retp(cfactor + 0.5*vec(abinv)'vec(sigma*abinv)-ln(abs(det(abinv))));
  endif;
endp;

proc lkgrad(theta);
  local binvt,temp,g1,g2;
  if z2==0;
    a=reshape(z1*theta+w1,k,k)';
    retp(z1'vec(sigma*a-inv(a')));
  elseif z1==0;
    b=reshape(z2*theta+w2,k,k)';
    binvt=inv(b)';
    retp(z2'vec(binvt-binvt*sigma*moment(binvt,0)));
  else;
    a=reshape(z1*theta+w1,k,k)';
    b=reshape(z2*theta+w2,k,k)';
    binvt=inv(b)';
    temp=sigma*a*moment(binvt,0);
    g1=z1'vec(temp-inv(a'));
    g2=z2'vec(binvt-binvt*a'temp);
    retp(g1+g2);
  endif;
endp;


proc lkinfo(theta);
  local ind,ainv,binvt,bibit,bainv,i11,i21,i22;
  ind=reshape(seqa(1,1,k)+seqa(0,k,k)',k*k,1);
  if z2==0;
    a=reshape(z1*theta+w1,k,k)';
    ainv=inv(a);
   retp(z1'(submat(ainv'.*.ainv,ind,0)+(eye(k).*.moment(ainv,0)))*z1);
  elseif z1==0;
    b=reshape(z2*theta+w2,k,k)';
    binvt=inv(b)';
    if scalmiss(binvt);
      retp(binvt);
    else;
      bibit=moment(binvt,0);
      retp(z2'((bibit.*.eye(k))+submat(binvt.*.binvt',ind,0))*z2);
    endif;
  else;
    a=reshape(z1*theta+w1,k,k)';
    b=reshape(z2*theta+w2,k,k)';
    ainv=inv(a);
    binvt=inv(b)';
    bibit=moment(binvt,0);
    bainv=b*ainv;
    i11=z1'(submat(ainv'.*.ainv,ind,0)+(bibit.*.moment(bainv,0)))*z1;
    i21=z2'((bibit.*.bainv)+submat(binvt.*.ainv,ind,0))*z1;
    i22=z2'((bibit.*.eye(k))+submat(binvt.*.binvt',ind,0))*z2;
    retp(i11-i21'-i21+i22);
  endif;
endp;

