[R] Non-Linear Regression (Cobb-Douglas and C.E.S)

[Sundar Dorai-Raj]

WilDscOp wrote:
> Dear all,
> 
> For estimating Cobb-Douglad production Function [ Y = ALPHA * 
> (L^(BETA1)) * (K^(BETA2))  ], i want to use nls function (without 
> linearizing it). But how can i get initial values?
> 
> ------------------------------------
>  > options(prompt=" R> " )
>  R> Y <- c(59.6, 63.9, 73.5, 75.6, 77.3, 82.8, 83.6, 84.9, 90.3, 80.5, 
> 73.5, 60.3, 58.2, 64.4, 75.4, 85, 92.7, 85.4, 92.3, 101.2, 113.3, 107.8, 
> 105.2, 107.1, 108.8, 131.4, 130.9, 134.7, 129.1, 147.8, 152.1, 154.3, 
> 159.9) # production
>  R> L <- c(39.4, 41.4, 43.9, 43.3, 44.5, 45.8, 45.9, 46.4, 47.6, 45.5, 
> 42.6, 39.3, 39.6, 42.7, 44.2, 47.1, 48.2, 46.4, 47.8, 49.6, 54.1, 59.1, 
> 64.9, 66, 64.4, 58.9, 59.3, 60.2, 58.7, 60, 63.8, 64.9, 66) # employment
>  R> K <- c(236.2, 240.2, 248.9, 254.5, 264.1, 273.9, 282.6, 290.2, 
> 299.4, 303.3, 303.4, 297.1, 290.1, 285.4, 287.8, 292.1, 300.3, 301.4, 
> 305.6, 313.3, 327.4, 339, 347.1, 353.5, 354.1, 359.4, 359.3, 365.2, 
> 363.2, 373.7, 386, 396.5, 408) # capital
>  R> klein <- cbind(Y,L,K)
>  R> klein.data<-data.frame(klein)
>  R> coef(lm(log(Y)~log(L)+log(K)))
> # i used these linearized model's estimated parameters as initial values
> (Intercept)      log(L)      log(K)
>  -3.6529493   1.0376775   0.7187662
>  R> nls(Y~ALPHA * (L^(BETA1)) * (K^(BETA2)), data=klein.data, start = 
> c(ALPHA=-3.6529493,BETA1=1.0376775,BETA2=0.7187662), trace = T)
> 6852786785 :  -3.6529493  1.0376775  0.7187662
> 1515217 :  -0.02903916  1.04258165  0.71279051
> 467521.8 :  -0.02987718  1.67381193 -0.05609925
> 346945.7 :   -0.5570735  10.2050667 -10.2087997
> Error in numericDeriv(form[[3]], names(ind), env) :
>         Missing value or an Infinity produced when evaluating the model
> ------------------------------------
> 
> 1. What went wrong? I think the initial values are not good enough: How 
> can i make a grid search?
> 

I think you meant this:

nls(Y~ALPHA * (L^(BETA1)) * (K^(BETA2)), data=klein.data, start = 
c(ALPHA=exp(-3.6529493),BETA1=1.0376775,BETA2=0.7187662), trace = TRUE)

Note that you neglected to exponentiate ALPHA.


> 2. How can i estimate C.E.S Production Function [  Y = GAMA * 
> ((DELTA*K^(-BETA)) + ((1-DELTA)*L^(-BETA)))^(-PHI/BETA)  ] using the 
> same data? How to get the initial value?
> 

This one is tougher. I set initial values DELTA = 0.5 and BETA = -1 and 
solved the log(Y) problem.

# log(Y) = log(GAMMA) + PHI * log(0.5 * K + 0.5 * L)
DELTA <- 0.5
BETA <- -1.0
start <- coef(lm(log(Y) ~ log(I(0.5 * K + 0.5 * L))))
GAMMA <- exp(start[1])
PHI <- -BETA * start[2]

However, I still couldn't get nls to even take one step due to singular 
gradient at the initial values. So I tried optim instead

ces.opt <- function(par, Y, L, K) {
   delta <- par[1]
   beta <- par[2]
   gamma <- par[3]
   phi <- par[4]
   Y0 <- delta * K^(-beta) + (1 - delta) * L^(-beta)
   Yhat <- gamma * Y0^(-phi/beta)
   sum((Y - Yhat)^2)
}

start <- optim(c(DELTA, BETA, GAMMA, PHI), ces.opt,
                Y = Y, L = L, K = K,
                method = "BFGS",
                control = list(maxit = 1000))
Y0hat <- with(start, par[1] * K^(-par[2]) + (1 - par[1]) * L^(-par[2]))
Yhat <- with(start, par[3] * Y0hat^(-par[4]/par[2]))
plot(Yhat ~ Y)

The fit wasn't to shabby. So I used the final values from the output 
from optim to feed to nls.

nls(Y ~ GAMMA * (DELTA * K^(-BETA) + (1 - DELTA) * 
L^(-BETA))^(-PHI/BETA), data=klein.data, start = c(DELTA = start$par[1], 
BETA = start$par[2], GAMMA = start$par[3], PHI = start$par[4]), trace = 
TRUE)

This still didn't work. Perhaps somebody else can explain this, though.


HTH,

--sundar