load("E:/3/dur_john.Rdata")
del=c(13,14,16,19,36,44,45,50,51,56,59,62,66,68,69,72,78,81,82,91,111,114,118)
data=data[-del,]#ɸѡһȥ

data=data[data$NONOIL==1,] #ɸѡ

sk=log(data$IONY)
n=log(data$POPGRO+5)#nʵԭn+g+delta
y=log(data$GDP85)
co1=sk-n
co2=co1^2

#basic Solow-CES model
#ԭеķ̣4
#unrestricted
library(sandwich)
library(lmtest)
lm1=lm(y~sk+n+co2)
summary(lm1)
coeftest(lm1,vcov=vcovHC(lm1,type="HC0"))

#restricted
#һʹDELTAֶϵȽ׼
lm2=lm(y~co1+co2)
summary(lm2)
b2=coeftest(lm2,vcov=vcovHC(lm2,type="HC0"))
#b2Ϊ󣬵һΪϵֵڶΪȽ׼
b2
#alphasigmaԵֵ
alpha=b2[2,1]/(1+b2[2,1])
alpha
sigma=1/(1-b2[3,1]*(1-alpha)^2*2/alpha)
sigma
#ֱӵnlsֱӵõϵȽ׼ȻҪֶ
lm2=lm(y~co1+co2)
lm3=nls(y~constant+alpha/(1-alpha)*co1+1/2*(sigma-1)/sigma*alpha/(1-alpha)^2*co2,start=list(constant=5,alpha=0.3,sigma=1.5),trace=F)
summary(lm3)

#׼
#Աx
X=cbind(1,co1,co2)
#ŸȾcc
cc=matrix(c(0,0,0,0),nrow=2,ncol=2)
cc[1,1]=1/(1-alpha)^2
cc[2,1]=(sigma-1)*(1-alpha^2)/(2*sigma*(1-alpha)^4)
cc[2,2]=alpha/(2*(1-alpha)^2*sigma^2)
#ɲвЭomegaݹʽ2õolsgammaЭm
omega=diag(lm2$resid^2)
m=solve(t(X)%*%X)%*%t(X)%*%omega%*%X%*%solve(t(X)%*%X)
#ݹʽ1õthetaЭvarԽϵֵŵõӦ׼
var=solve(cc)%*%m[-1,-1]%*%solve(t(cc))
se=diag(var)^0.5
alpha_se=se[1]
alpha_se
sigma_se=se[2]
sigma_se
