R语言之非线性回归xt9.5

第9章 非线性回归

9.5 表9-13数据中GDP和投资额K都是用定基居民消费价格指数CPI缩减后的值，1978年的价格指数为100。（C-D生产函数y=AKαLβ）
（1）用线性化的乘性误差项模型拟合C-D生产函数。
（2）用非线性的最小二乘拟合加性误差项的C-D生产函数。
（3）对线性化回归检验自相关，如果存在自相关则用自回归方法改进。
（4）对线性化回归检验共线性，如果存在共线性则用岭回归方法改进。

年份 t CPI GDP k l
1978 1 100.00 3624.1 1377.9 40152
1979 2 101.90 3962.9 1446.7 41024
1980 3 109.54 4124.2 1451.5 42361
1981 4 112.28 4330.6 1408.1 43725
1982 5 114.53 4623.1 1536.9 45295
1983 6 116.82 5080.2 1716.4 46436
1984 7 119.97 5977.3 2057.7 48197
1985 8 131.13 6836.3 2582.2 49873
1986 9 139.65 7305.4 2754.0 51282
1987 10 149.85 7983.2 2884.3 52783
1988 11 178.02 8385.9 3086.8 54334
1989 12 210.06 8049.7 2901.5 55329
1990 13 216.57 8564.3 2975.4 64749
1991 14 223.94 9653.5 3356.8 65491
1992 15 238.27 11179.9 4044.2 66152
1993 16 273.29 12673.0 5487.9 66808
1994 17 339.16 13786.9 5679.0 67455
1995 18 397.15 14724.3 6012.00 68065
1996 19 430.12 15782.8 6246.50 68950
1997 20 442.16 16840.6 6436.00 69820
1998 21 438.62 17861.6 6736.10 70637
1999 22 432.48 18975.9 7098.90 71394
2000 23 434.21 20604.7 7510.50 72085
2001 24 437.25 22256.0 8567.30 73025
2002 25 433.75 24247.0 9764.90 73740

rm(list=ls())

# ---- xt9.5 数据中GDP和投资额K都是用定基居民消费价格指数CPI缩减后的值，1978年的价格指数为100 ----
# ---- （1）用线性化的乘性误差项模型拟合C-D生产函数。 ----
data9.5 <- read.csv('D:/rwork/应用回归/习题数据/表9-13.csv',head=TRUE)
t <- data9.5[,2]
y <- data9.5[,4]
ly <- log(data9.5[,4])
lk <- log(data9.5[,5])
ll <- log(data9.5[,6])
model1 <- lm(ly~lk+ll,data9.5)
summary(model1)
anova(model1)
# 线性化模型：A=0.17，α=0.80，β=0.40。  #A=exp(-1.78545)



# ---- （2）用非线性的最小二乘拟合加性误差项的C-D生产函数。 ----
model2 <- nls(y~A*((k^a)*(l^b)),data9.5,start=list(A=2,a=0.9,b=0.3),
              lower=c(0,0,0),upper=c(10000,100,100),algorithm='port',
              control=nls.control(maxiter=1000,tol=1e-1000))
summary(model2)
# 非线性化模型：A=0.41，α=0.87，β=0.27。



# ---- （3）对线性化回归检验自相关，如果存在自相关则用自回归方法改进。 ----
## 检验自相关-DW检验 ----
library(lmtest)
dwtest(model1,alternative='two.sided') #DW检验
# DW=0.72，存在自相关。
# 可知DW值为0.72，P值=3.1e-05，查DW表，n=25,k=2,显著性水平α=0.05，
#  得dL=1.29,dU=1.45,由于DW=0.72<dL=1.20,知DW值落入正相关区域，即残差序列存在正的自相关。



## 用迭代法处理序列相关，并建立方程 ----
# 一次迭代ly(t)'=ly(t)-ρ*ly(t-1)，lk(t)'=lk(t)-ρ*lk(t-1)，ll(t)'=ll(t)-ρ*ll(t-1)
# 自相关系数ρ^=1-DW/2=1-0.72/2=0.64，计算的ly',lk',ll'
rho_hat <- 1-0.72/2 #自相关系数ρ^
n <- length(ly)
lyy <- ly[2:n]-rho_hat*ly[1:n-1]
lkk <- lk[2:n]-rho_hat*lk[1:n-1]
lll <- ll[2:n]-rho_hat*ll[1:n-1]
model3 <- lm(lyy~lkk+lll)
summary(model3) #回归分析
anova(model3) #方差分析表
# 用迭代法得到A=0.40，α=0.73，β=0.53。  #A=exp(-0.90814)



# ---- （4）对线性化回归检验共线性，如果存在共线性则用岭回归方法改进。 ----
## 检验共线性-方差扩发因子法 ----
library(car)
vif(model1)
#  所有自变量对应的VIF全部大于10，所以自变量之间存在共线性。



## 检验共线性-特征根判定法 ----
XX <- cor(data.frame(ly,lk,ll)) #计算样本相关阵
kappa(XX,exact=TRUE) #exact=TRUE时，精确计算条件数，否则近似计算条件数。
# 根据条件数100<=k=870.1625<=1000，说明自变量之间存在较强的多重共线性。



## 建立ly对lk、ll的岭回归。 ----
datas <- data.frame(scale(data.frame(ly,lk,ll))) 
#对样本数据进行标准化处理并转换为数据框的格式存储
library(MASS)
ridge9.5 <- lm.ridge(ly~lk+ll,data=datas,lambda=seq(0,8,0.2))
# 做岭回归，对于标准化后的数据模型不包含截距项，其中lambda为岭参数k的所有取值

beta <- coef(ridge9.5) #将所有不同岭参数所对应的回归系数的结果赋给beta
beta
# 结果中，第1列为岭参数k，其取值范围为0到3，步长0.1，共有31个k值。
# 第2至3列是数据标准化后的岭回归系数，其中第1行k=0的数值就是普通最小二乘估计的标准化回归系数。

## 绘制岭迹图
k <- ridge9.5$lambda #将所有岭参数赋给k
plot(k,k,type='n',xlab='岭参数k',ylab='岭回归系数',ylim=c(0.4,1.15))
# 创建没有任何点和线的图形区域
linetype <- c(1:2)
char <- c(18:19)
for (i in 1:2)lines(k,beta[,i+1],type='o',lty=linetype[i],pch=char[i],cex=0.75) #画岭迹线
legend(4.5,1,inset=0.5,legend=c('lk','ll'),cex=0.8,pch=char,lty=linetype) #添加图例
# legend(locator(i),inset=0.5,legend=c('lk','ll'),cex=0.8,pch=char,lty=linetype) #添加图例 # 注意运行时鼠标点击双击图片

# 两个自变量的方差扩大因子皆大于13,且条件数大于870,存在严重多重共线性。取岭回归参数=5，得岭估计A=0.00083, α=0.48, β=1.13。

参考课本：应用回归分析（R语言版），何晓群编著