/*
**
**      LSTAR model of the Australian unemployment rate
**
*/

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

format /rd 10,3;




/**     Read in quarterly data for Australia from March 1971 to December 2008       **/



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

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


unemployment   = xlsreadm("ausu.xls", "b1:b152", 1,vls);            /**     unemployment rate            **/


/**     Generate LSTAR variables      **/

dy = trimr(unemployment,1,0) - trimr(unemployment,0,1);             /**     Dependent variable: dy(t)       **/

y_lag = trimr(unemployment,0,1);                                    /**     Explanatory variable: y(t-1)    **/

dy_lag1 = trimr(dy,3,1);                                          /**     Explanatory variable: dy(t-1)   **/
dy_lag2 = trimr(dy,2,2);                                          /**     Explanatory variable: dy(t-2)   **/
dy_lag3 = trimr(dy,1,3);                                          /**     Explanatory variable: dy(t-3)   **/
dy_lag4 = trimr(dy,0,4);                                          /**     Explanatory variable: dy(t-4)   **/


dy = trimr(dy,4,0);         /**     Readjust these variables as a result of dulag   **/
y_lag = trimr(y_lag,4,0);


t = rows(dY);



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

proc lnlt(b,y);

     local u,mue1,mue2,alpha1,alpha2,gam,c,g,sig2,lnl;

     mue1 = b[1];                                              /**     Define the parameters               **/
     alpha1 = b[2];


     mue2 = b[3];
     alpha2  = b[4];

     gam = b[5];
     c      = b[6];

     sig2  = b[7]^2;


    g      = 1./(1 + exp(-gam*(dy_lag1 - c)));

     u     = dy - (mue1 + alpha1*y_lag + b[8]*dy_lag1  + b[9]*dy_lag2  + b[10]*dy_lag3 + b[11]*dy_lag4) 
                - (mue2 + alpha2*y_lag + b[12]*dy_lag1 + b[13]*dy_lag2 + b[14]*dy_lag3 + b[15]*dy_lag4).*g;

     lnl = - 0.5*ln(2*pi) - 0.5*ln(sig2) - 0.5*u.^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
**/

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


theta_0 = 0.18 | -0.03 | 0.2 | -0.5 | 2 | -0.1 | 0.12 | 1*ones(8,1);


theta_0 = {
0.1839941074978994 
-0.0291767302348609 
-0.2764566105735068 
0.1168435291737226 
0.1506415237143992 
0.3381271551483930 
0.2595982566281941 
-0.2293473702070491 
0.2796543465420748 
-0.2470307094856789 
-0.7799268183024534 
-0.8337906047214504 
7.1707315254410080 
0.5975982331610694 
0.2811750295052042 


};



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

    __title = "LSTAR Model";

    { theta,a,g,vcov,retcode } = maxprt(maxlik(dy,0,&lnlt,theta_0));


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)      **/



/**     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,dy) );

endp;


mt = gradp(&lnlt1,theta);

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;



