/*   This is file QuaidEst.gau.  It estimates the semiflexible QUAIDS model with the Slutsky 
	substitution matrix (evaluated at the sample mean) restricted to be symmetric negative 
	semi-definite and of rank k(<=n-1), using the Cholesky decomposition method. 
	"NSystem" determines which system is estimated.
1. First, estimate by OLS the NSystem*(n-1)-parameter Unrestricted QUAIDS(0) ["Unrestricted" 
	in that lambda's are raw rather than deflated by beta(p); k=0, so all gamma's=0].
	(This is also ML, so LF value is calculated directly.)
2. Next, use OLS params. as starting values for ML estimation of symmetric concave 
	QUAIDS, system for k=0, with the exact alpha(p) price index, for k=0.  
	(Estimation proceeds by iterating on beta(p).)
3. Repeat step 2. for selected values of k between 1 and n-1.  (At every stage, use the 
	parameter values from the immediately preceding stage as starting values.)
N.B. StQUAIDS(n-1) with symmetry but without NSD is not calculated

(c) J.P. Neary 1999

*/
new; /*************************************************************************/
/*          USER-SPECIFIED VARIABLES: ECONOMIC DATA AND PARAMETERS  */
n=11;                       /* n = # goods        N.B. This and next two overridden if TESTIF=1 */
nyear=1980;
kklast=n;  /* kklast=1+klast, where klast is max. value of k estimated: 
      k cannot exceed n-1, so kklast cannot exceed n */
 /*************************************************************************/
/*    USER-SPECIFIED VARIABLES WHICH CONTROL OPERATION OF PROGRAM  */
library optmum, gearylib, pgraph;
NSystem=2;     /* NSystem = 1 for HAIDS; =2 for AIDS; =3 for QUAIDS */
ntol=10^(-14);	/* ntol:  Tolerance level for convergence of ML algorithm; =10^(-14) in a full run */
NExact=1;        /* NExact=0 if Stone version is to be estimated; 1 for Exact version */
TESTIF=0;        /* TESTIF=1 on a test run (using small data set only)  */
printif=0; 	        /* printif=1 for printing of output from SORTDATA */
printOLS=1;     /* printOLS= 1 for printing of OLS Summary Table */
printif9=1;        /* printif9=1 for printing of Summary Table of each ML iteration*/ 
GRAPHIF=0;    /* graphif=0: no graph; >1: print graphs */
SCREEN ON;
/*************************************************************************/
IF NSystem == 1;
    save path = c:\me\research\geary2\empirics\gaussout\HAIDS;
    screen on; print  "Estimates of HAIDS System";
    goto READNOW;
 ELSEIF NSystem == 2;
    save path = c:\me\research\geary2\empirics\gaussout\AIDS;
    screen on; print  "Estimates of AIDS System";
    goto READNOW;
ENDIF;
save path = c:\me\research\geary2\empirics\gaussout\QUAIDS;
    screen on; print  "Estimates of QUAIDS System";
READNOW:
/*************************************************************************/
{ m,p,q } = readdata(n,nyear,testif);  /* Call PROC READDATA to read raw data and print
   heading:  m=#countries; p and q are n-by-m matrices of prices and quantities */
countrynum=seqa(1,1,m);
/**************************************************************************/
/*  INITIALISE MATRICES     */
alphaOLS=zeros(n,1);  /* n-by-1 vector of OLS alpha's:  #n to be estimated residually */
betaOLS=zeros(n,1);               /* n-by-1 vector of OLS beta's:  #n to be estimated residually */
lamOLS=zeros(n,1);  /* n-by-1 vector of OLS lambdas:  #n to be estimated residually */
teeOLS=zeros(n,3);               /* n-by-3 vector of OLS t-values:  #n a dummy (always zero) */
rsqOLS=zeros(n,1);               /* n-by-1 vector of OLS R-squareds:  #n a dummy (always zero) */
resid=zeros(n-1,m);    /* (n-1)-by-m matrix of OLS residuals */
alpha=zeros(n-1,1);    /* (n-1)-by-1 vector of ML alpha's:  n'th to be estimated residually */
beta=zeros(n-1,1);      /* (n-1)-by-1 vector of ML beta's:  n'th to be estimated residually */
lambda=zeros(n-1,1);  /* (n-1)-by-1 vector of ML lambdas: n'th to be estimated residually */
ggamma=zeros(n-1,n-1);   /* (n-1)-by-(n-1) matrix of ML gamma's:  n'th to be estimated residually */
tau=zeros(n-1,n-1);            /* (n-1)-by-(n-1) matrix of Cholesky decomposition of Slutsky matrix */
thetalst=zeros((n-1)*(n+6)/2,1);     /* vector of parameter values from previous loop:
	big enough to hold full theta vector for maximum value of k; 
	normally to be reinitialised at end of ML loop) */
LFML=zeros(n,1);   /* n-by-1 vector of maximised log likelihood values for each value of k */
Timerun=zeros(n,1);   /* n-by-1 vector of time taken in each round of ML DO-loop */
Nparm=zeros(n,1);         /* n-by-1 vector of no. of params. corresponding to each value of k */
alphaALL=zeros(n,n-1);    /* n-by-(n-1) matrix of ML alpha's for different values of k */
betaALL=zeros(n,n-1);      /* n-by-(n-1) matrix of ML beta's for different values of k */
lamALL=zeros(n,n-1);      /* n-by-(n-1) matrix of ML lambdas for different values of k */
gammaALL = zeros(n*(n-1),n-1);  /* matrix of (n-1)-by(n-1) gammas for all n values of k */
/***********************************************************************/
/* Call PROC SORTDATA to sort p and q by z, scale p and z by sample mean 
     and calculate: pprimeq, z, w, w_trim, poverz and poz_trim */
{p, q, pprimeq, z, zraw, w, w_trim, poverz, poz_trim}=SORTDATA(p,q,printif);
logprn=log(p[1:n-1,.])-log(p[n,.]);
pStone = sumc(w.* log(p)); /* Compute log of Stone price index: m-by-1 vector, 
	giving column sums of n-by-m matrix of element-by-element product of w and log(p) */
logofyS = log(z) - pStone;  /* Log of nominal expenditure deflated by Stone price index */
/***********************************************************************/
/*  FIRST STEP:  CALCULATE OLS ESTIMATES OF k=0 SPECIFICATION */
aa1=100;  /* Note:  Residuals scaled by aa1, to avoid over-small digits and tiny det.  30.11.98 */
IF NSYSTEM==1;  /* HAIDS: OLS ests. of alpha are just the average budget shares */
   alphaOLS= meanc(w'); 
   resid = aa1*(w[1:n-1,.] - alphaOLS[1:n-1,1]*ones(1,m)); /* Calculate (n-1)-by-m matrix of resids. */
   GOTO CALCLF;
 ELSEIF NSystem == 2;  /* AIDS:  Calculate OLS ests. of unrestricted Stone-AIDS[0] */
   _olsres=1;  /* Instruct proc OLS to compute OLS residuals */
   i=1; 
  DO while i<=n-1; 
   SCREEN Off; print "OLS ests. of 2-param. Unrestricted Stone-AIDS eqtn. with k=0, for good" i;
   {vname,mom, b, stb,vc,stderr,sigma,cx,rsq,residi,dwstat}=OLS(0,(w[i,.])',log(z)-pStone);
   alphaOLS[i,1]=b[1,1];     betaOLS[i,1]=b[2,1];  
   teeOLS[i,1:2]=(b./stderr)';   /* teeOLS[i,3] a dummy in AIDS case */    rsqOLS[i,1]=rsq;
   resid[i,.]=aa1*residi'; /* read residuals from this eqtn. into matrix of residuals */
    i=i+1; 
   ENDO; 
   GOTO ENDQUAID;
 ELSEIF NSystem == 3; /*  QUAIDS:  Calculate OLS ests. of unrestricted Stone-QUAIDS[0]  */
   _olsres=1;  /* Instruct proc OLS to compute OLS residuals */
   x=(log(z)-pStone)~ ((log(z)-pStone)^2);  /* LHS variables for each regression  */
   /* x=(log(z)-pStone)~ ((log(z)-pStone)^2)./(10^pStone);  /* LHS variables: Stone deflation case */  */
  i=1; 
  DO while i<=n-1; 
     SCREEN Off; print "OLS ests. of 3-param. Unrestricted Stone-QUAIDS eqtn. with k=0, for good" i;
     {vname,mom, b, stb,vc,stderr,sigma,cx,rsq,residi,dwstat}=OLS(0,(w[i,.])',x);
     alphaOLS[i,1]=b[1,1];     betaOLS[i,1]=b[2,1];      lamOLS[i,1]=b[3,1];
     teeOLS[i,.]=(b./stderr)';   rsqOLS[i,1]=rsq;
     resid[i,.]=aa1*residi'; /* read residuals from this eqtn. into matrix of residuals */
     i=i+1; 
  ENDO; 
  lamOLS[n,1]=-ones(n-1,1)'*lamOLS[1:n-1,1];   /* calculate lambda for good #n by subtraction */
ENDIF;
ENDQUAID:  /* Next two tidying-up calculations work for both AIDS and QUAIDS cases */
alphaOLS[n,1]=1-ones(n-1,1)'*alphaOLS[1:n-1,1];   /* calculate alpha for good #n by subtraction */
betaOLS[n,1]=-ones(n-1,1)'*betaOLS[1:n-1,1];   /* calculate beta for good #n by subtraction */
/***********************************************************************/
CALCLF:   /* Calculate Maximised Value of OLS Likelihood Function Directly  */  
LFvalue=  -0.5*m*(n-1)*(log(2*pi)+1) - 0.5*m* log (det(resid*resid'/m)) + m*(n-1)*log(aa1); 
	 /* Value of likelihood function at maximum (corrected for scaling of residuals)*/
/***********************************************************************/
/*                      PRINT OLS SUMMARY INFORMATION               */
IF PRINTOLS==1;
SCREEN On; print "Unrestricted OLS Estimates with k=0";
format /rd 9,4;  print "   group #     a-OLS    t-value     b-OLS    t-value lambda-OLS t-value R squared" 
    seqa(1,1,n)~alphaOLS~teeOLS[.,1]~betaOLS~teeOLS[.,2]~lamOLS~teeOLS[.,3]~rsqOLS;  
SCREEN On;  print "LF value: " LFvalue "; (scaling factor for residuals =" aa1 ")";
ENDIF;
/******************************************************************************/
/******************************************************************************/
/*  MAIN LOOP: CALCULATE ML ESTIMATES FOR EACH VALUE OF  k  */
/* Loop over kk = 1+k, where k is the number of non-zero rows in tau; 
        hence rank of tau and Slutsky matrix; maximum of k: n-1; so max of kk: n */
kk=1; 
ntaulast=0;  /* Value of ntau in "previous" round: see B. and C. below */
/******************************************************************************/
DO WHILE kk<=kklast; /* MASSIVE DO LOOP!!!  */
/******************************************************************************/
/*  SET CONTROL PARAMETERS FOR THIS ROUND */
Timenow=hsec;
k=kk-1;
nfirst=NSystem*(n-1);   /* number of parameters other than tau: n-1 each of alpha, beta, 
	lambda [Only alpha in HAIDS; only alpha and beta in AIDS ] */
ntau=0.5*k*(2*n-k-1);  /* number of non-zero elements in tau matrix */
Nparm[kk,1]=nfirst+ntau; /* Note: index of Nparm (=# parameters) is kk, NOT k */
/******************************************************************************/
/* CREATE STARTING VALUES FOR DEMAND PARAMETERS IN THIS ROUND */
/* Vector of starting values for this round: theta0 is (nfirst+ntau)-by-1, as follows:
    A. nfirst [=NSystem*(n-1)] values comprising: (n-1) each of alpha, beta and lambda:
              OLS values when k=0, otherwise ML values from previous round.
    B. (n-k) [=ntau-ntaulast] NEW elements of tau; arbitrary, but do NOT set equal to zero initially!
    C. ntaulast values of tau from previous round (so, of course, irrelevant for k=0 or 1) 
B. and C. combined: Use Nparm[kk-1,1] parameters from previous round       */
/******************************************************************************/
IF k==0;   /* IFF k=0: (a) OLS starting values; (b) ntau=0 so theta0 contains no tau terms */ 
    beta=betaOLS[1:n-1,1];  /* OLS starting values for beta */
    IF NSystem == 1;     /*  HAIDS:  OLS starting values for alpha only  */
          theta0= alphaOLS[1:n-1,1] ;    
          goto MLSTART;
    ELSEIF NSystem == 2;   /* AIDS:  OLS starting values for alpha and beta only */
          theta0= alphaOLS[1:n-1,1] | betaOLS[1:n-1,1] ;   /* k=0: [A] only */
          goto MLSTART;
    ELSEIF NSystem == 3;   /* NSystem = 3,  QUAIDS:  OLS starting values for alpha, beta and gamma */
          theta0= alphaOLS[1:n-1,1] | betaOLS[1:n-1,1] | lamOLS[1:n-1,1];   /* k=0: [A] only */
          goto MLSTART;
    ENDIF;
ENDIF;

theta0=thetalst[1:Nparm[kk-1,1],1] | 0.1+zeros(n-k,1);   /* k>0: [A], [B] & [C] */
/***********************************************************************/
/* SUB-LOOP 1: CALCULATE ML ESTS. FOR CURRENT VALUE OF k */
MLSTART:  /***********************************************************/
screen on; format /rd 3,0; print "Next Round:  kk  k  ntau  ntaulast  rows(theta0): "
        kk  k  ntau  ntaulast  rows(theta0);  format /rd 5,3; 
print  "theta0': " theta0';  
/***********************************************************************/
aa=100;
betaofp=10^((logprn)'*beta);   /* m-by-1 vector of beta(p) functions */
logofy = logofyS;  /* Initial value of m-by-1 vector of log expenditure deflated by alpha(p) 
	(needed in PROC LFSQUAk) uses Stone price index as deflator */
/***********************************************************************/
/* CALCULATE ML ESTS. OF EXACT QUAIDS FOR CURRENT VALUE OF k */
/*  N.B.  No need to reset betaofp: just start with values from last Stone round.
   This means that first run through next loop is a dummy, since parameter values are
    unchanged.  (Done this way to avoid having to recalculate alpha_n and gammabit twice.) */
kkk=1;
LFlast=0;  
DO WHILE kkk<=100;  /* Iterate over beta(p) and logofy up to 100 times.   */
   SCREEN Off;  {theta,f,g,retcode}=OPTPRT(OPTMUM(&LFSQUAk,theta0)); 
   LFML[kk,1]=  -0.5*m*(n-1)*(log(2*pi)+1) - 0.5*m* log (f) + m*(n-1)*log(aa); 
   screen on; format /rd 8,4; print kkk LFML[kk,1];
   IF (LFML[kk,1]-LFLast)^2<= ntol;    /* Test for convergence of Likelihood Function */
	diff=theta-theta0;
	IF diff'*diff <= ntol;	/* Test for convergence of whole theta vector */
		GOTO ENDExact;
	ENDIF;
   ENDIF;
   LFLast=LFML[kk,1];
   betaofp=10^((logprn)'*beta);   /* update vector of beta(p) functions */
   /************************************************************************/
   /*  Update vector of alpha(p) functions (See PROC LOGALFAP) and deflate expenditure by it */
   logofy = log(z) - logalfap(p);  /* Value of m-by-1 vector of log deflated 
	expenditure (needed in PROC LFSQUAk) uses exact deflator [alpha(p)]  */
   /************************************************************************/
   theta0=theta;      /* Current values of theta are starting values for next kkk round */
   kkk=kkk+1;
ENDO;
ENDExact:
/************************************************************************/
WRAPUP:     /* Finish off the current round */
eigs=eigh(-tau'*tau);  /* Calculate eigenvalues of Slutsky substitution matrix at sample mean*/
/*  Store parameters for current value of k in summary matrices */
alphaALL[kk,.]=alpha';    
betaALL[kk,.]=beta';     
lamALL[kk,.]=lambda';
nfirstga = (kk-1)*(n-1) + 1;
nlastga = kk*(n-1);
gammaALL[nfirstga:nlastga,.] = ggamma;
Timerun[kk,1]=(hsec-Timenow)/6000;
nprint: /***********************************************************************/
/*            PRINT FINAL SUMMARY INFORMATION FOR CURRENT VALUE OF k               */
IF PRINTIF9==1;
SCREEN Off; format /rd 2,0; print "MaxLik Estimates of Semiflexible QUAIDS of Rank " k "with " nfirst+ntau "Parameters";
SCREEN On; format /rd 6,4; print "tau: " tau; 
print "ML alpha, beta, lambda and gamma: " alpha~beta~lambda~ggamma; 
SCREEN On;  format /rd 9,6;  print "LF value: " LFML[kk,1] "; (scaling factor for residuals =" aa ")";
SCREEN On; format /rd 7,5; print "Eigenvalues of Slutsky Substitution Matrix:\L" eigs'; 
ENDIF;
/************************************************************************/
ENDLOOP:
ntaulast=ntau;  /* Store ntau from this round, to determine dimension of theta0 in NEXT round */
thetalst=theta;	 /* Store parameter values from this round for use in theta0 in NEXT round */
kk=kk+1;    /* Update index of DO-loop */
ENDO;   /* END OF MASSIVE DO-LOOP!  */  
/************************************************************************/
/*  Save parameters to disk, for subsequent use by program GaiaAIDS.GAU */
save alphaALL;   save betaALL;  save lamALL;  save gammaALL;
/************************************************************************/
/*   PRINT OVERALL SUMMARY  */
format /rd 2,0; print "Summary for all " kk-1 "Specifications Estimated \L     k     # params.  LF value   Time Taken";
format /rd 9,4; print seqa(0,1,n)~Nparm~LFML~Timerun;
SCREEN On; print "k:    alpha's:"; format /rd 7,4; print seqa(0,1,n)~alphaALL;
print "k:       beta's:"; format /rd 7,4; print seqa(0,1,n)~betaALL;
print "k:       lambda's:"; format /rd 7,4; print seqa(0,1,n)~lamALL;
/***********************************************************************/
/*   PRINT GRAPH OF VALUES OF LOG LIKELIHOOD */
graphset; screen on;   /*  SET UP GRAPH CONTROLS */
title("Maximised Log Likelihood for Exact QUAIDS");
xtics(0,90,10,0);   ytics(430,490,10,0);   /*  Set axis tick marks   */
_plctrl = {1};  /* Draw line and plot symbol for every point on the LFML curve */
xy(Nparm,LFML);
/************************************************************************/
finish:
end;


/*  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
PROC LFSQUAk:  Calculates Likelihood Function for semiflexible QUAIDS of rank k */
PROC LFSQUAk(theta);  /*  Input: theta is an (nfirst+ntau)-by-1 vector; Returns:  ssqresid  */
local kk, nend, nstart, wijhat, resid;  
/* Globals used:  n, m, k, NSystem, nfirst, alpha, beta, lambda, ggamma, tau, betaofp, logofy,
	logprn, aa, w_trim */     /* Reminder:  nfirst = NSystem*(n-1)  */
/************************************************************************/
/*  Read nfirst parameters form theta into alpha, beta and lambda vectors; defaults for them are zero */
alpha=theta[1:n-1,1];  /* (n-1)-by-1:  alpha vector equals first n-1 entries in theta */ 
IF NSystem >=2;
    beta=theta[n:2*n-2,1];  /* (n-1)-by-1; NSystem>1 only:  beta equals next n-1 entries in theta */
    IF NSystem == 3;
        lambda=theta[2*n-1:3*n-3,1];  /* (n-1)-by-1; NSystem=3 only:  lambda vector equals next n-1 entries in theta */
    ENDIF;
ENDIF;
/************************************************************************/
IF k==0;   /* When (and only when) k=0, ntau=0 so theta contains no tau terms */ 
	GOTO OMEGAHAT;
ENDIF;
kk=1;     /* Start loop to construct tau from remaining (and last) ntau elements of theta */
DO while kk<=k;    /* Read elements of theta into row kk of tau matrix */
	nend = nfirst+0.5*kk*(2*n-kk-1);  	/* Index of element of theta corresponding to 
		last entry in row # kk of tau */
	nstart=nend-(n-kk-1);  			/* Index of element of theta corresponding to 
		first entry in row # kk of tau */
	tau[kk,kk:n-1]=theta[nstart:nend,1]';
	kk=kk+1;
ENDO;
ggamma = - tau'*tau - alpha*alpha' + diagrv(zeros(n-1,n-1),alpha); /* (n-1)-by-(n-1):  gamma 
	matrix formed from tau and alpha */ 
OMEGAHAT:  /*****************************************************************/
/*  Calculate wijhat */
/*   Deflate lambda's by betaofp (i.e., the exact deflator) */
wijhat=alpha*(ones(m,1))' + beta*(logofy)' + lambda*(((logofy)^2)./betaofp)'
      + ggamma*logprn;   /*  (n-1)-by-m matrix of predicted budget shares */     
/************************************************************************/
resid=aa* (w_trim-wijhat) ;  /* calculate (n-1)-by-m matrix of residuals;
              Note:  Residuals scaled by aa, to avoid over-small digits and tiny det.  30.11.98 */
retp(det(resid*resid'/m));
endp;