ECONOMICS 951
NUMERICAL METHODS FOR
ECONOMICS AND ECONOMETRICS

Assignment 2.0

Write a function that estimates the parameters of the selection problem

(1)     y* = x1 beta1 + u1
(2)     m* = x2 beta2 + u2
(3)     m  = (m*>0)
(4)     y  = y* if m and y = .  if !m  ('.'=missing)
(5)     (u1,u2) jointly normal; stdev(u1)=sigma; var(u2) = 1; corr(u1,u2)=rho;
(6)     theta = beta1 
                beta2 
                sigma > 0
                rho, |rho|<1

(1)     y* = b11 + b12x2 + b13x3 + u1
(2)     m* = b21 + b22x2 + b23x3 + b24x4 + u2
That is, x4 is excluded from (1) a priori.  

Details

1. In your main function:

    
       Read the data using loaddta()
       Split the input matrix into y,m, and X 
       Concatenate a column of 1s in X
       Create a vector s which describes the structural system of equations:  
           s[j]=0 means X[.,j] appears in neither x1 nor x2;
           s[j]=1 means X[.,j] appears in x1 only;
           s[j]=2 means X[.,j] appears in x2 only;
           s[j]=3 means X[.,j] appears in both x1 and x2; 
       call hecklee defined below
       

2. Write these functions (you have to fix syntax and fill in details):

   
       probit(const beta2, const adFunc, const avScore, const amHessian) {
           // Returns the probit likelihood function for equation (3).  
           // Supply the data through global variables.
           decl pb;
           pb = probn(-X2*beta2);
           return (pb).^m (1-pb).^(1-m) ???;
           } 
       
       imill(const d) {
           //d a column vector 
           //Return a column vector of corresponding inverse Mill's ratios.  
           return -densn(d) ./ (1-probn(d)) ???;
           }
           
       vfiml(const theta) {
           // compute the full information likelihood function for the model (1)-(6)
           // return N x 1 vector containing the likelihood for each observation.  
           // Supply the data through global variables.          
           pb = probn();
           dn = ???;
           return pb .* dn ???;
           }
   
        numfd(const func,const theta,const gmatrix) {
            //   theta is k x 1
            //   func(theta) is Nx1 
            //   gmatrix on output contains the is N x k matrix of 
            //     numerical first derivatives of func. 
            //   get gradient step length from a global variable (eps2?)
            return ;
            }
   
       fiml(const theta,const like, const Score, const BHHH) {
               //    return likelihood, score and hessian for the model (1)-(6)  
           vlk = vfiml(const theta);
           like=sumc(vlk);
           if (Score!=0) or  (BHHH!=0) { 
               numfd(vfiml,theta,gmat);
               if Score!=0 { &Score=sumc(gmat) };
               if BHHH!=0 { &BHHH = gmat' * gmat};
               };
           return 1 ???;        
           }
   
       hecklee(const fiml, const y,const m, const X,const s,const theta, const vhtheta) {
       //    If  fiml=0 return the Heckman two-step estimates of theta. 
       //    If  fiml=1 then return the full-information maximium likelihood estimates of theta
       //               and the BHHH estimates of its variance matrix.
       //    (Of course estimates of the variance should be returned in the two-step casee, 
       //        but this is left as optional.)
           process the input vectors
           get starting values for beta2. (use OLSC?)
           MaxBFGS(probit, beta2, lkmax, 0, 0)
           lamhat = imilla(-X2*beta2)
           olsc(y, X1~lamhat, beta1star, xxinv)
           &theta = stuff from beta2, beta1star, and xxinv. 
           if fiml=1 { 
               use theta from the 2-step method as starting values to maximize fiml().
               use some combination of MaxSimplex(), MaxBFGS(), and MaxNewton(). 
               &BHHH = BHHH estimate of the variance matrix.
               }
           return 1(?)
           }