/*
**
**      Example taken from Greene, Table 13.1, 3rd edition.
**
**      Test an AR(1) model of investment using the conditional MLE.
**
*/

library maxlik;
#include maxlik.ext;
maxset;
gausset;

cls;

format /rd 14,8;

output file=invest.out reset;


/**     Read in quarterly data for the US over the period 1957q1 to 2007q2, Source: IFS  **/


chdir c:\diskette\econbook\answers\computer;

vls = reshape(error(0),9,1);

r10yr  = xlsreadm("investment","aj2:ic2",1,"")';        /**     government bond 10-year (% pa)  **/
invest = xlsreadm("investment","aj3:ic3",1,"")';        /**     nominal investment (sa)         **/
gdp    = xlsreadm("investment","aj4:ic4",1,"")';        /**     nominal gdpt (sa)               **/
r3yr   = xlsreadm("investment","aj5:ic5",1,"")';        /**     government bond 3-year (% pa)   **/
cpi    = xlsreadm("investment","aj6:ic6",1,"")';        /**     consumer price index)           **/
tbill  = xlsreadm("investment","aj7:ic7",1,"")';        /**     Treasury bill rate (% pa)       **/


r10yr~invest~gdp~r3yr~cpi~tbill;


/**     Generate variables: data start in 1958q1 as a result of loosing 4 observations form computing annual inflation rate **/

lri = ln(trimr(invest./cpi,1,0));                       /**     real investment (logs)          **/

lry = ln(trimr(gdp./cpi,1,0));                          /**     real income (logs)              **/

inf = ln(trimr(cpi,1,0)./trimr(cpi,0,1));           /**     inflation (% pa)                **/

rint = trimr(r10yr/100,1,0) - 4*inf;                          /**     real interest rate  (% pa)      **/



t = rows(lri);


y = lri;
x = ones(t,1)~lry~rint;


itoggle = 0.0;      /**     toggle to ensure that rho is in the unit circle: choose 0 for estimation and 1 for standard errors  **/


/**     Define the unconstrained conditional log of the likelihood at each observation            **/

proc lnlt1(b,y);

     local u,v,beta0,beta1,beta2,rho1,sig2,lnl,lnl_1;

     beta0 = b[1];                                              /**     Define the parameters               **/
     beta1 = b[2];
     beta2 = b[3];

     if itoggle == 0.0; rho1  = tanh(b[4]); else; rho1 = b[4]; endif;   /**     Restriction to ensure that estimate of rho is in unit circle    **/

     sig2  = abs(b[5]);

     u     = lri - beta0 - beta1*lry - beta2*rint;

     v     = trimr(u,1,0) - rho1*trimr(u,0,1);

     lnl   = - 0.5*ln(2*pi) - 0.5*ln(sig2) - 0.5*v.^2/sig2;     /**     Log-likelihood for t=2,3,...        **/

     retp( lnl );

endp;



/**     Define the constrained log of the likelihood at each observation            **/

proc lnlt0(b,y);

     local u,v,beta0,beta1,beta2,rho1,sig2,lnl,lnl_1;

     beta0 = b[1];                                              /**     Define the parameters               **/
     beta1 = b[2];
     beta2 = b[3];
     rho1  = 0.0;
     sig2  = abs(b[4]);

     u     = lri - beta0 - beta1*lry - beta2*rint;

     v     = trimr(u,1,0) - rho1*trimr(u,0,1);

     lnl   = - 0.5*ln(2*pi) - 0.5*ln(sig2) - 0.5*v.^2/sig2;     /**     Log-likelihood for t=2,3,...        **/

     retp( lnl );

endp;


/**     Algorithm choices:
                            steep
                            bfgs
                            dfp
                            newton
                            bhhh
                            prcg
**/

/**     Covariance choices:
                            nocov
                            info
                            xprod
                            hetcon
**/

bols = y/x;

e  = y - x*bols;
s2 = meanc(e.^2);

rols = trimr(e,1,0)/trimr(e,0,1); 


_max_options = { bfgs info }; /* Newton-Rap algorithm, Hessian std errors */
                                /* Leave a space between brackets and names */


/**     Estimate the unconstrained model using Newton-Raphson and compute Hessian se  **/

__title = "Unconstrained Model";

theta_0 = bols | rols | s2 ;    /**     Construct starting estimates based on OLS       **/

{ theta,a1,g1,vcov1,retcode } = maxprt(maxlik(y,0,&lnlt1,theta_0));


itoggle = 1; theta_0 = theta[1:3] | tanh(theta[4]) | theta[5];               //  Reestiamte to derive correct standard errors

{ theta1,a1,g1,vcov1,retcode } = maxprt(maxlik(y,0,&lnlt1,theta_0));

lnl1 = (t-1)*a1;                            /* Unconstrained log-likelihood        */



/**     Estimate the constrained model using Newton-Raphson and compute Hessian se  **/

__title = "Constrained Model";

theta_0 = bols | s2 ;

{ theta0,a0,g0,vcov0,retcode } = maxprt(maxlik(y,0,&lnlt0,theta_0));

lnl0 = (t-1)*a0;                           /*  Constrained log-likelihood           */


/**     Likelihood ratio test       **/

lr = -2*(lnl0 - lnl1);


print "";
print "Log of the likelihood (unconstrained) = " lnl1;
print "Log of the likelihood (constrained)   = " lnl0;
print "";
print "Likelihood ratio test                 = " lr;
print "p-value                               = " cdfchic(lr,1);
print "";




/**     Wald test       **/


r = 0 ~ 0 ~ 0 ~ 1 ~ 0;

q = 0;

w = (r*theta1 - q)'inv(r*vcov1*r')*(r*theta1 - q);

print "";
print "vcov (unconstrained)" vcov1;
print "";

print "Wald test                      = " w;
print "p-value                        = " cdfchic(w,1);
print "";



/**     Lagrange Multiplier test (based on numerical opg matix)              **/

proc lnlt(b);                               /**     This procedure is used to compute the gradients numerically of the unconstrained model      **/

    retp( lnlt1(b,y) );

endp;


theta0 = theta0[1 2 3] | 0.0 | theta0[4];              /**     Fill in constrained parameter estimates     **/

gmat = gradp(&lnlt,theta0);

g = sumc(gmat);

j = gmat'gmat;

lm = g'*inv(j)*g;

print "";
print "gradient evaluated at contrained estimates" g';
print "";
print "numerical opg" j;
print "";
print "var-covar    " inv(j);
print "";
print "LM test (based on numerical opg)      = " lm;
print "p-value                               = " cdfchic(lm,1);
print "";  



/**     Lagrange Multiplier test (based on Gauss-Newton version)            **/

y = ri;

x = ry ~ rint;

b = y/x;                /**     Regression equation under the null          **/

e = y - x*b;            /**     Residuals under the null                    **/

z = ones(t-1,1) ~ trimr(x,1,0) ~ trimr(e,0,1);

y = trimr(e,1,0);       /**     Redefine residuals                          **/

b = y/z;                /**     Second stage regression estimates           **/

e = y - z*b;            /**     OLS residuals from second stage regression  **/

r2= 1 - e'e/(y-meanc(y))'(y-meanc(y));              /**  r-squared          **/

lm_gn  = (t-1)*r2;

print "";
print "LM test (based on Gauss-Newton)       = " lm_gn;
print "p-value                               = " cdfchic(lm_gn,1);
print "";  



/**     Lagrange Multiplier test (based on first order autocorrelation of residuals)            **/

y = ri;

x = ry ~ rint;

b = y/x;                /**     Regression equation under the null          **/

e = y - x*b;            /**     Residuals under the null                    **/

r1 = acf(e,1,0);        /**     First order autocorrelatin of residuals     **/

lm_r1  = (t-1)*r1^2;

print "";
print "LM test (based on first order ar1)    = " lm_r1;
print "p-value                               = " cdfchic(lm_r1,1);
print "";  


output off;



/*****************************

Log of the likelihood (unconstrained) =   568.60435951 
Log of the likelihood (constrained)   =   302.95099289 

Likelihood ratio test                 =   531.30673323 
p-value                               =     0.00000000 


vcov (unconstrained)
    0.60240102    -0.06164089     0.00593086    -0.00247712     0.00000001 
   -0.06164089     0.00986039    -0.00082931     0.00008730     0.00000000 
    0.00593086    -0.00082931     0.00270554    -0.00001505     0.00000000 
   -0.00247712     0.00008730    -0.00001505     0.00002020     0.00000000 
    0.00000001     0.00000000     0.00000000     0.00000000     0.00000000 

Wald test                      = 48892.08236665 
p-value                        =     0.00000000 


gradient evaluated at contrained estimates   -0.00134516     0.00308028    -0.00000087   191.49811958     0.00042411 

numerical opg
71743.45904537 313276.37477493  2147.77601315 -1178.07421889 -203589.45160164 
313276.37477493 1371968.27276910  9481.53316869 -5568.21633156 -965085.72225473 
 2147.77601315  9481.53316869   135.62226970   -91.10252410 -16179.12806237 
-1178.07421889 -5568.21633156   -91.10252410   598.08613056 74230.60458520 
-203589.45160164 -965085.72225473 -16179.12806237 74230.60458520 13752845.29121373 

var-covar    
    0.00532796    -0.00122413     0.00042318    -0.00008153    -0.00000609 
   -0.00122413     0.00028267    -0.00019850     0.00001937     0.00000138 
    0.00042318    -0.00019850     0.01586430     0.00011206     0.00001039 
   -0.00008153     0.00001937     0.00011206     0.00506980    -0.00002708 
   -0.00000609     0.00000138     0.00001039    -0.00002708     0.00000024 

LM test (based on numerical opg)      =   185.91736761 
p-value                               =     0.00000000 


LM test (based on Gauss-Newton)       =   187.80287390 
p-value                               =     0.00000000 


LM test (based on first order ar1)    =   187.69536129 
p-value                               =     0.00000000 


**********************************/

