/* Simfin5.prg
** GAUSS ECSUR algorithm using simulated data
** Parke Wilde pwilde@econ.ag.gov
** 3-13-99
** Modified 9-23-99
*/

/*****************************************************************************
CAUTION:
Please use this code feely, but entirely at your own risk.
I would be eager to hear any suggestions or corrections, and also to hear how this code has been used.
This code replaces an earlier program ECSUR that was posted on the Gauss archives and withdrawn due to an error (an interesting error, perhaps: it turns out that there is no known iterative SUR method for estimating this model with "unbalanced" panel data.
My apologies to one known user of the previous program, and any unknown users!)
The current program uses direct maximum likelihood estimation through the MAXLIK routines.
E-mail pwilde@econ.ag.gov.

OVERVIEW:
Estimates a system of M equations that are linear in the parameters.
Error terms have an error-components structure, with one component that is shared by all members t=1,...,Tn of household n (or all time periods observed for individual n), and one component that is idiosyncratic for each member.
There are n=1,...,N households.
Both error components may be correlated across equations, so there are two cross-equation variance-covariance matrices to estimate, each of which is rather like the one estimated in more traditional SUR models.
Currently, this code assumes the same variables appear in each equation.
Unlike the usual SUR, this case does not produce estimates that are equivalent to OLS.


REFERENCES:
Baltagi (1995), Econometric Analysis of Panel Data, Wiley and Sons.
Magnus (1982), "Multivariate Error Components Analysis of Linear and Nonlinear Regression Models by Maximum Likelihood," Journal of Econometrics 19:239-285.
Wilde, McNamara, and Ranney (1999), "The Effects of Income and Food Programs on Dietary Quality: A Seemingly Unrelate Regression Analysis with Error Components," American Journal of Agricultural Economics 81(4):201-213.
Pollak and Wales, Demand System Specification and Estimation, Oxford Press.  Pp. 134-137.

CONTENTS:
Section 1   -- Makes simulated data
Section 2   -- Defines Globals used later
Section 3   -- Procedures for likelihood functions and gradients
Section 4   -- Calls estimation procedures and prints simulation results

PROGRAM STRUCTURE:
Section 4 calls maximum likelihood estimation routine MAXLIK.  MAXLIK is instructed to use likelihood functions and gradients in section 3: callsur and calldsur, respectively.  Because the data may be unbalanced (having a different number of members in each hh), these routines make use of subroutines liksur and dliksur, which calculate the likelihood function and gradients for "balanced" subsets of the data where each hh has the same number of members.

The global vector panindex keeps track of how many hh in the data matrix dat have each possible size: row i of panindex reports how many hh in the data have i members.  In a real application, the user should sort the data matrix by hh size (so that the first block of data contains all hh with one member, the next block all hh with 2 members, and so forth).  The simulated data contain this structure already.  Zeros are permitted in panindex (there need not be positive numbers of hh with each number of members).  In a real application, the user must provide panindex (and proof read it, there is no error trapping provided for errors in panindex).  The reason for this somewhat severe data structure is to limit the number of times that sub-routines liksur and dliksur must be called on each iteration of the maximization algorithm.

Global vector yind identifies which columns of dat contain the dependent variables.  Global vector xind identifies which columns of dat contain the independent variables (including a column of ones if an intercept is needed).  Global scalar tind identifies which column of dat contains the hh size.  Each row of dat contains data for one member of one hh.  If that member is not the last member of the hh, the next row will contain the next member of the same hh.  For example, if panindex={2,0,2} and yind={1,2} and xind={3,4} and tind=5, the data matrix dat has the structure:
hh#:  member#:     dat:
1     1            yvar.eq1  yvar.eq2  1  xvar  1
2     1            yvar.eq1  yvar.eq2  1  xvar  1
3     1            yvar.eq1  yvar.eq2  1  xvar  3
3     2            yvar.eq1  yvar.eq2  1  xvar  3
3     3            yvar.eq1  yvar.eq2  1  xvar  3
4     1            yvar.eq1  yvar.eq2  1  xvar  3
4     2            yvar.eq1  yvar.eq2  1  xvar  3
4     3            yvar.eq1  yvar.eq2  1  xvar  3

Parameter vector b contains all parameters, including regression parameters beta and the unique elements of the two variance-covariance matrices.  Several globals identify the elements of parameter vector b.  Vector parbeta reports which elements of b contain beta, in order: first 1,...,k for equation 1, then 1,...,k for equation 2, and so forth until equation m.  Vector parv reports which elements of b contain the unique elements of the variance-covariance matrix for the individual error component.  Vector parm reports which elements of b contain the unique elements of the variance-covariance matrix for the hh error component.  For example, in the simulation there are two equations with an intercept and two other variables each, plus the variance-covariance matrices.  So, parameter vector b has 12 elements: parbeta={1,2,3,4,5,6}, parv={7,8,9}, and parm={10,11,12}.

If user has installed MAXLIK (including maxlik.ext), this program is set up
to run the simulation.
*****************************************************************************/

/***********************************  SECTION 1  ****************************/
new;
library maxlik;
#include d:\pwgauss\maxlik\src\maxlik.ext;

/* Panel Index */
/* This vector is needed because data may be unbalanced. */
/* Each row i of this vector tells how many hh have exactly i members. */
/* The simulation makes data with 1000 hh with 2 individuals, and 1000 hh */
/* with 4 individuals, but those parameters may be adjusted. */
panindex={0, 1000, 0, 1000};

/*
proc newsim:  Makes a block of simulated data, where all hh have same # indivs
Inputs:
  b     K by M matrix of "true" parameters
  sigv  M by M matrix of "true" variance-covariance for indiv error component
  sigm  M by M matrix of "true" variance-covariance for hh error component
  t     number of individuals in each hh
  n     number of hh with t individuals
Outputs:
  dat   data matrix with nt observations
  yind  pointer to the column of dat containing dependent variable y
  xind  pointer to the columns of dat containing independent variables x
  tind  pointer to the columns of dat containing tcol, # of ind in each hh */
proc (4) = newsim(b,sigv,sigm,t,n);
  local k,m,tobs,x,mu,v,y,tcol,dat,yind,xind,tind;
  tobs=t*n;
  k=rows(b)-1;
  m=cols(b);
  if (m/=rows(sigv))+(m/=rows(sigm))+(m/=cols(sigv))+(m/=cols(sigm))>0;
    print "Warning.  Invalid input to proc simecsur.";
  endif;
  x=ones( tobs ,1)~5*rndu( tobs ,k);     /* regressors */
  if maxc(maxc(sigm)) >0;
    mu=rndn(n,m)*chol(sigm);          /* error components "mu" */
  endif;
  mu=mu.*.ones(t,1);
  print "dim mu";
  print rows(mu) cols(mu);
  v=rndn(tobs,m)*chol(sigv);         /* error components "v" */
  print "dim v";
  print rows(v) cols(v);
  y=x*b+mu+v;                                  /* y variables */
  print "dim y";
  print rows(y) cols(y);
  tcol=t*ones(tobs,1);
  dat=x~y~tcol;                                /* full data set "dat" */
  yind=seqa(k+2,1,m);  /* indices identifying y variables */
  xind=seqa(1,1,k+1);  /* indices identifying x variables */
  tind=k+m+2;          /* index identifying size of cross-section */
retp(dat,yind,xind,tind);
endp;

/* Simulated "True" Parameters */
/* True v-cv matrix for error component mu shared by all members of hh */
truesigm={
2    0.5,
0.5    2  };
/* True v-cv matrix for individual error component v */
truesigv={
2    0.25,
0.25   2  };
/* True regression parameter matrix */
trueb={2 1,4 2,6 4};

/* Create data set "dat" with 6000 rows (1000*2 rows + 1000*4 rows) */
dat2={};
dat4={};
{dat2,yind,xind,tind}=newsim(trueb,truesigv,truesigm,2,panindex[2]/2);
{dat4,yind,xind,tind}=newsim(trueb,truesigv,truesigm,4,panindex[4]/4);
dat=dat2|dat4;

/***********************************  SECTION 2  ****************************/
/* Globals for sample size, etc. */
tobs=rows(dat);
m=2;
k=2;
mv=(m*(m+1)/2);

/* Make duplication matrix used in proc dliksur below */
/* A computational convenience explained in Magnus appendix */
s=seqa(1,1,mv);
proc xpand(s);
  retp(vec(xpnd(s)));
endp;
dup=round(gradp(&xpand,s));
print "dimension dup";
print rows(dup) cols(dup);

/* Globals saying which elements of parameter vector are beta, which are v-cv matrices */
parbeta=seqa(1,1,(k+1)*m);
parv=seqa((k+1)*m+1,1,mv);
parm=seqa((k+1)*m+mv+1,1,mv);

/* Starting vector - Uses untrue parameters beta to demonstrate convergence */
untruth={0.1,0.1,0.2,-0.1,-0.1,-0.2};
bstart=(vec(trueb)+untruth)|vech(truesigv)|vech(truesigm);

/***********************************  SECTION 3  ****************************/
/* Likelihood proc for all hh with Tn members */
/* This proc must be called separately, first for all hh with 1 member, then for
** all hh with 2 members, and so forth.  It is called by proc callsur below */
proc liksur(b,data);
  local t, tobs, n, beta, sigv, sigmu, sig1, jbar, e, utmp,
        ncount ,rownum,invsigv, invsig1, const, llik, newdetsig, llik2;
  t=data[1,tind];
  tobs=rows(data);
  n=tobs/t;
  beta=reshape(b[parbeta],m,k+1)';

  sigv= xpnd( b[parv] );
  sigmu=xpnd( b[parm]);
  sig1=sigv+t*sigmu;
  invsigv=invpd(sigv);
  invsig1=invpd(sig1);

  jbar=ones(t,t)/t;
  e=eye(t)-jbar;

  const= t*m*ln(2*pi) + ln(det(sig1))+(t-1)*ln(det(sigv));

  llik={};
  ncount=0;
  do until ncount==n;
    rownum=seqa((t*ncount+1),1,t);
    utmp=data[rownum,yind]-data[rownum,xind]*beta;
    llik=llik|diag(utmp*invsig1*utmp'*jbar+utmp*invsigv*utmp'*e);
    ncount=ncount+1;
  endo;
  llik=-0.5*(const/t+llik);
  retp(llik);
endp;

/* Gradient */
proc dliksur(b,data);
  local t, tobs, n, beta, sigv, sigmu, sig1, jbar, e, utmp,
        ncount ,rownum,invsigv, invsig1, invsig, dbet , xtmp , sig1til, sigvtil,
        xbig, ubig, db, dall /* , dup */;
  t=data[1,tind];
  tobs=rows(data);
  n=tobs/t;
  beta=reshape(b[parbeta],m,k+1)';
  /* dup={1 0 0, 0 1 0, 0 1 0, 0 0 1};/*"duplication matrix" see Magnus appendix*/ */
  sigv= xpnd( b[parv] );
  sigmu=xpnd( b[parm]);
  sig1=sigv+t*sigmu;
  invsigv=invpd(sigv);
  invsig1=invpd(sig1);

  jbar=ones(t,t)/t;
  e=eye(t)-jbar;

  invsig=invsig1.*.jbar + invsigv.*.e;

  dall={};
  ncount=0;
  do until ncount==n;
    /* rownum=seqa((t*ncount+1),1,t);
    xtmp=data[rownum,xind];
    utmp=data[rownum,yind]-xtmp*beta; */
    xtmp=data[(t*ncount+1):(t*ncount+t),xind];
    utmp=data[(t*ncount+1):(t*ncount+t),yind]-xtmp*beta;
    xbig=eye(m).*.xtmp;
    ubig=vec(utmp);

    dbet=(xbig'*invsig*ubig)';

    sig1til=invsig1*(utmp'*jbar*utmp-sig1)*invsig1;
    sigvtil=invsigv*(utmp'*e   *utmp-(t-1)*sigv)*invsigv;

    db=dbet~((1/2)*dup'*vec(sig1til+sigvtil))'~((t/2)*dup'*vec(sig1til))';
    /* dall=dall|(1/t)*ones(t,1)*db; */
    dall=dall|db;
    ncount=ncount+1;
  endo;
  retp(dall);
endp;

/* Calls likelihood proc liksur.  Is called by maxlik. */
proc callsur(b,data);
  local tc, llik, first, last;
  llik={};
  tc=0;
  do until tc==rows(panindex);
    tc=tc+1;
    if panindex[tc]>0;
      last=sumc(panindex[1:tc]);
      first=last-panindex[tc]+1;
      /* print "sur first" first;
      print "sur last" last; */
      llik=llik|liksur(b,data[first:last,.]);
    endif;
  endo;
retp(llik);
endp;

/* Calls gradient proc dliksur.  Is called by maxlik. */
proc calldsur(b,data);
  local tc, dlik, first, last;
  dlik={};
  tc=0;
  do until tc==rows(panindex);
    tc=tc+1;
    if panindex[tc]>0;
      last=sumc(panindex[1:tc]);
      first=last-panindex[tc]+1;
      /* print "dsur first" first;
      print "dsur last" last; */
      dlik=dlik|dliksur(b,data[first:last,.]);
    endif;
  endo;
retp(dlik);
endp;

/***********************************  SECTION 4  ****************************/
/*  Estimate Using Maximum Likelihood */
  _max_Algorithm = 4;
  _max_GradCheckTol = 1e-2;
_max_maxIters=30;
/*  _max_covPar=2;  */
_MAX_GRADPROC=&calldsur;
{betamax,f,g,cov,retcode}= maxlik(dat,0,&callsur,bstart);
call maxprt(betamax,f,g,cov,retcode);

/* Print Results for Comparison to True Parameters */
print "betamax";
print reshape(betamax[parbeta],m,k+1)';
print "trueb";
print trueb;
print "sigv";
print xpnd(betamax[parv]);
print "truesigv";
print truesigv;
print "sigm";
print xpnd(betamax[parm]);
print "truesigm";
print truesigm;

