/*
**
**      Estimate an AR(1) model of investment by exact and conditional maximum likelihood
**
*/

library maxlik,pgraph;
#include maxlik.ext;
maxset;
gausset;
graphset;

cls;

format 10,6;

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;



/**
        Choose the algorithm        

                algor = 0 for exact MLE

                algor = 1 for conditional MLE

**/


algor = 0.0;

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


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

proc lnlt(b,y);

     local u,v,beta0,beta1,beta2,rho1,sig2,lnl,lnl_0,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 = - 0.5*ln(2*pi) - 0.5*ln(sig2) + 0.5*ln(1 - rho1^2) - 0.5*(u[1] - 0).^2/(sig2/(1 - rho1^2));    /**     Log-likelihood for t=1              **/

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

     if algor == 0.0; 

        lnl = lnl_0 | lnl_1;

        else;

        lnl = lnl_1;

    endif;

    retp( lnl );

endp;


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

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

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


bols = y/x;

e  = y - x*bols;

s2 = meanc(e.^2);

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

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



if algor == 0; 

    /**     Estimate the model by exact MLE using Newton-Raphson and compute Hessian se  **/

    __title = "Exact MLE of the AR(1) Model";

    { b,a,g,vcov,retcode } = maxprt(maxlik(y,0,&lnlt,bstart));

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

    { b,a,g,vcov,retcode } = maxprt(maxlik(y,0,&lnlt,bstart));

    else;

    /**     Estimate the model by conditional MLE using Newton-Raphson and compute Hessian se  **/

    __title = "Conditional MLE of the AR(1) Model";

    { b,a,g,vcov,retcode } = maxprt(maxlik(y,0,&lnlt,bstart));

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

    { b,a,g,vcov,retcode } = maxprt(maxlik(y,0,&lnlt,bstart));

endif;

h = -inv(vcov);                     /**     Compute the hessian from the output of MAXLIK (make sure that the INFO option is used to compute the standard errors in MAXLIK)      **/


print "";
print "Average log-likelihood =" a;

print "";
print "Log-likelihood         =" a*(t-1);

print "";
print "Variance-covariance matrix" vcov;

print "";
print "Hessian" h;




/**     Compute covariances (White)   **/


proc lnlt1(b);                      /**     This is a dummy procedure that is needed by gradp to compute the numerical gradients of the log-likelihood     **/

    retp( lnlt(b,y) );

endp;


mt = gradp(&lnlt1,b);

j = mt'mt;                          /**     Outer product of gradients              **/

qmle_white = inv(h)*j*inv(h);



/**     Compute covariances (Newey-West)   **/

t = rows(mt);

lmax = floor( 4*(t/100)^(2/9) );

        w  = mt'mt;

        tau=1;

        do until tau>lmax;
           wtau = mt[(tau+1):t,.]'mt[1:(t-tau),.];
           w    = w + (1.0-tau/(lmax+1))*(wtau + wtau');
           tau  = tau + 1;
        endo;

qmle_nw = inv(h)*w*inv(h);


print "";
print "Standard errors (Hessian)    " sqrt(diag(inv(-h)))';
print "Standard errors (OPG)        " sqrt(diag(inv(j)))';
print "Standard errors (White)      " sqrt(diag(qmle_white))';
print "Standard errors (Newey-West) " sqrt(diag(qmle_nw))';
print "      based on lags equal to " lmax;

output off;



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

===============================================================================
                         Exact MLE of the AR(1) Model                          
===============================================================================
 MAXLIK Version 5.0.6                                      3/15/2011   6:02 pm
===============================================================================

return code =    0
normal convergence

Mean log-likelihood         2.81479
Number of cases     201

Covariance matrix of the parameters computed by the following method:
Inverse of computed Hessian

Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
------------------------------------------------------------------
P01             -2.8908        0.4105   -7.042   0.0000      0.0000
P02              1.3041        0.0878   14.858   0.0000      0.0000
P03             -0.0475        0.0531   -0.894   0.3711      0.0000
P04              3.3855        0.7149    4.736   0.0000      0.0000
P05              0.0002        0.0000    9.994   0.0000      0.0000

Correlation matrix of the parameters
   1.000  -0.885   0.122  -0.368   0.026
  -0.885   1.000  -0.137   0.435  -0.031
   0.122  -0.137   1.000  -0.062   0.004
  -0.368   0.435  -0.062   1.000  -0.071
   0.026  -0.031   0.004  -0.071   1.000

Number of iterations    39
Minutes to convergence     0.00230

===============================================================================
                         Exact MLE of the AR(1) Model                          
===============================================================================
 MAXLIK Version 5.0.6                                      3/15/2011   6:02 pm
===============================================================================

return code =    0
normal convergence

Mean log-likelihood         2.81479
Number of cases     201

Covariance matrix of the parameters computed by the following method:
Inverse of computed Hessian

Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
------------------------------------------------------------------
P01             -2.8908        0.4103   -7.045   0.0000      0.0000
P02              1.3041        0.0877   14.866   0.0000      0.0000
P03             -0.0475        0.0531   -0.894   0.3711      0.0000
P04              0.9977        0.0033  305.842   0.0000      0.0000
P05              0.0002        0.0000    9.994   0.0000      0.0000

Correlation matrix of the parameters
   1.000  -0.885   0.122  -0.367   0.026
  -0.885   1.000  -0.137   0.434  -0.031
   0.122  -0.137   1.000  -0.062   0.004
  -0.367   0.434  -0.062   1.000  -0.071
   0.026  -0.031   0.004  -0.071   1.000

Number of iterations    0
Minutes to convergence     0.00000

Average log-likelihood =  2.814794 

Log-likelihood         =562.958882 

Variance-covariance matrix
  0.168365  -0.031862   0.002665  -0.000491   0.000000 
 -0.031862   0.007695  -0.000639   0.000124   0.000000 
  0.002665  -0.000639   0.002824  -0.000011   0.000000 
 -0.000491   0.000124  -0.000011   0.000011   0.000000 
  0.000000   0.000000   0.000000   0.000000   0.000000 

Hessian
-27.494026 -115.047840   0.189335  75.596206  -1.215280 
-115.047840 -644.024812 -28.859113 2185.469170  12.620038 
  0.189335 -28.859113 -360.849413 -20.938720  -0.407641 
 75.596206 2185.469170 -20.938720 -116502.292525 -1064817.791839 
 -1.215280  12.620038  -0.407641 -1064817.791839 -2397059946.251238 

Standard errors (Hessian)      0.410323   0.087721   0.053145   0.003262   0.000020 
Standard errors (OPG)          0.331681   0.064692   0.055737   0.007586   0.000020 
Standard errors (White)        0.558734   0.123986   0.051609   0.002386   0.000023 
Standard errors (Newey-West)   0.541312   0.121821   0.053196   0.002394   0.000023 
      based on lags equal to   4.000000 





===============================================================================
                      Conditional MLE of the AR(1) Model                       
===============================================================================
 MAXLIK Version 5.0.6                                      3/15/2011   6:01 pm
===============================================================================

return code =    0
normal convergence

Mean log-likelihood         2.84302
Number of cases     201

Covariance matrix of the parameters computed by the following method:
Inverse of computed Hessian

Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
------------------------------------------------------------------
P01             -4.4265        0.7805   -5.671   0.0000      0.0000
P02              1.5134        0.0994   15.229   0.0000      0.0000
P03             -0.0638        0.0520   -1.226   0.2203      0.0000
P04              2.8860        0.3662    7.882   0.0000      0.0000
P05             -0.0002        0.0000   -9.994   0.0000      0.0000

Correlation matrix of the parameters
   1.000  -0.799   0.147  -0.714   0.000
  -0.799   1.000  -0.161   0.199   0.000
   0.147  -0.161   1.000  -0.065   0.000
  -0.714   0.199  -0.065   1.000   0.000
   0.000   0.000   0.000   0.000   1.000

Number of iterations    67
Minutes to convergence     0.00567

===============================================================================
                      Conditional MLE of the AR(1) Model                       
===============================================================================
 MAXLIK Version 5.0.6                                      3/15/2011   6:01 pm
===============================================================================

return code =    0
normal convergence

Mean log-likelihood         2.84302
Number of cases     201

Covariance matrix of the parameters computed by the following method:
Inverse of computed Hessian

Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
------------------------------------------------------------------
P01             -4.4265        0.7761   -5.703   0.0000      0.0000
P02              1.5134        0.0993   15.240   0.0000      0.0000
P03             -0.0638        0.0520   -1.226   0.2202      0.0000
P04              0.9938        0.0045  221.116   0.0000      0.0000
P05             -0.0002        0.0000   -9.994   0.0000      0.0000

Correlation matrix of the parameters
   1.000  -0.800   0.147  -0.710   0.001
  -0.800   1.000  -0.161   0.196   0.000
   0.147  -0.161   1.000  -0.064   0.000
  -0.710   0.196  -0.064   1.000  -0.001
   0.001   0.000   0.000  -0.001   1.000

Number of iterations    0
Minutes to convergence     0.00000

Average log-likelihood =  2.843022 

Log-likelihood         =568.604360 

Variance-covariance matrix
  0.602401  -0.061641   0.005931  -0.002477   0.000000 
 -0.061641   0.009860  -0.000829   0.000087   0.000000 
  0.005931  -0.000829   0.002706  -0.000015   0.000000 
 -0.002477   0.000087  -0.000015   0.000020   0.000000 
  0.000000   0.000000   0.000000   0.000000   0.000000 

Hessian
-39.912627 -214.400407  -0.297936 -3968.051770   2.700582 
-214.400407 -1259.622248 -32.212102 -20871.746787 -29.413621 
 -0.297936 -32.212102 -379.835586 -180.308169   0.740743 
-3968.051770 -20871.746787 -180.308169 -446031.915800 -6939.553298 
  2.700582 -29.413621   0.740743 -6939.553298 -2678758739.728256 

Standard errors (Hessian)      0.776145   0.099299   0.052015   0.004494   0.000019 
Standard errors (OPG)          0.685706   0.071633   0.053556   0.004798   0.000020 
Standard errors (White)        0.992235   0.140651   0.051790   0.004452   0.000020 
Standard errors (Newey-West)   1.025823   0.135230   0.056105   0.004386   0.000021 
      based on lags equal to   4.000000 

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

