/*  These data are taken from the SAS book Forecasting Examples, p. 3.
    The data is orginally from Kmenta, J. (1986), Elements of Econometrics
    second ed. New York: Macmillan Publishing Co., Inc.  This cross-
    sectional CES production function data set contains observations of
    the capital (K) and labor (L) and the outputs (Q) of 25 firms.
    This data is ideal for demonstrating the application of nonlinear
    least squares.  */

data ces;
   input id k l q @@;
   label k ='Capital Input'
         l ='Labor Input'
         q ='Output';
   datalines;
01 8 23 106.00  02 9 14  81.08  03 4 38  72.80  04 2 97  57.34
05 6 11  66.79  06 6 43  98.23  07 3 93  82.68  08 6 49  99.77
09 8 36 110.00  10 8 43 118.93  11 4 61  95.05  12 8 31 112.83
13 3 57  64.54  14 6 97 137.22  15 4 93  86.17  16 2 72  56.25
17 3 61  81.10  18 3 97  65.23  19 9 89 149.56  20 3 25  65.43
21 1 81  36.06  22 4 11  56.92  23 2 64  49.59  24 3 10  43.21
25 6 71 121.24
;

data ces;
  set ces;
  logq = log(q);

PROC NLIN DATA=ces best=8  method = gauss outest=test;
      
      PARMS A0 = 0 to 3 by 0.2
            alpha = 0 to 1 by 0.2
            h = 0.5 to 1.5 by 0.2
            rho = 0 to 0.9 by 0.1;

      MODEL logq = A0 - (h/rho)*log(alpha*k**(-rho) + (1-alpha)*l**(-rho)) ;

/* NOTE: If you choose not to provide SAS with the derivatives DER, but,
         do specify METHOD = GAUSS in the options of PROC NLIN, SAS will
         compute numerical derivatives and carry out the estimation. Results
         in that case may not necessarily be identical to the case when you
         do specify derivatives. See SAS/ETS Guide for further reference. */
     
  DER.A0 = 1;
  
  DER.alpha = -(h/rho)*(1/(alpha*k**(-rho) + (1-alpha)*l**(-rho)))*(k**(-rho) - l**(-rho));
  
  DER.h = -log(alpha*k**(-rho) + (1-alpha)*l**(-rho))/rho;
  
  DER.rho = log(alpha*k**(-rho) + (1-alpha)*l**(-rho))*(h/rho**2) +
                (h/rho)*(1/(alpha*k**(-rho) + (1-alpha)*l**(-rho)))*
                (alpha*k**(-rho)*log(k) + (1-alpha)*l**(-rho)*(log(l)));

RUN;

PROC NLIN DATA=ces outest=test method = gauss;
      
      PARMS A0 = 2 alpha = 0.5  h = 1 rho = 0.3;

      MODEL logq = A0 - (h/rho)*log(alpha*k**(-rho) + (1-alpha)*l**(-rho)) ;

  DER.A0 = 1;
  
  DER.alpha = -(h/rho)*(1/(alpha*k**(-rho) + (1-alpha)*l**(-rho)))*(k**(-rho) - l**(-rho));
  
  DER.h = -log(alpha*k**(-rho) + (1-alpha)*l**(-rho))/rho;
  
  DER.rho = log(alpha*k**(-rho) + (1-alpha)*l**(-rho))*(h/rho**2) +
                (h/rho)*(1/(alpha*k**(-rho) + (1-alpha)*l**(-rho)))*
                (alpha*k**(-rho)*log(k) + (1-alpha)*l**(-rho)*(log(l)));

RUN;

PROC NLIN DATA=ces outest=test method = gauss;
      
      PARMS A0 = 2 alpha = 0.5  h = 1 rho = 0.3;

      MODEL logq = A0 - (h/rho)*log(alpha*k**(-rho) + (1-alpha)*l**(-rho)) ;
      
RUN;

/* This part of the program is to print out the CORRELATION MATRIX of the
   Nonlinear Least Squares estimates.  */

proc print data=test;
run;

/*  Now convert the correlation matrix of the NLLS estimates into the 
    VARIANCE-COVARIANCE matrix of the NLLS estimates.  D is a diagonal
    matrix with the standard errors of the estimates on the diagonal.  
    CORR is the correlation matrix reported by the OUTEST option in
    Proc NLIN.  Therefore VCOV is the variance-covariance matrix of 
    the NLLS estimates.  */

proc iml;

D = {0.1120 0 0 0,
     0 0.1273 0 0,
     0 0 0.0501 0,
     0 0 0 0.2597};
             
CORR = {1 0.1743625 -0.6696781 -0.0652757,
        0.1743625 1 0.5813928 -0.9815590,
        -0.6696781 0.5813928 1 -0.6271445,
        -0.0652757 -0.9815590 -0.6271445 1};

VCOV = D*CORR*D;

/* Now let's set up the restriction matrix, R, for testing that the CES production
   function can be simplified to the Cobb-Douglas production function.  If h = 1 AND
   rho = 0 then the CES production function reduces to the Cobb-Douglas production
   function with constant returns to scale.  rho = 0 implies the Cobb-Douglas
   form but with increasing or decreasing returns to scale.   */ 
   
/* Proc IML code for computing WALD test statistic */

R = {0 0 1 0,
	 0 0 0 1};

q = {1,0};

theta = {2.4323,0.4064,0.8222,0.6043}; 

W = (R*theta - q)`*inv(R*VCOV*R`)*(R*theta - q);

critval = cinv(.95,2);
pval = 1 - probchi(W,2);

print "Results of Wald Test for Cobb-Douglas Production Function with CRS";

print "The Wald test statistic is :", W;

print "Chi-square (2) 5% critical value is:", critval;

print "P-value of Wald test is:", pval; 

quit;