引用信息：Gong B. Agricultural reforms and production in China: Changes in provincial production function and productivity in 1978–2015[J]. Journal of Development Economics, 2018, 132: 18-31.

在传统的效率分析中，通常使用随机前沿分析（SFA）和数据包络分析（DEA）进行测算效率。SFA 的生产函数通常设置为 Cobb-Dougles 或者 Translog 形式，并且引入随机误差项来解决外生冲击的影响，例如农业生产中的天气变化。DEA 则不假设生产前沿的具体函数形式，但没有到实际生产中具有随机性。在面板数据中，传统的 SFA 和 DEA 的劣势会被进一步放大。

有一些研究对 SFA 和 DEA 进行了改进，尝试解决面板数据中传统效率分析方法的问题，例如 TFE-SFA（True Fixed Effect SFA）、内生 SFA（Endog-SFA）等。然而，这些基于 SFA 方法的改进仍然未解决一个基础问题：不能考察随时间变化导致的投入产出关系的变化。 在农业领域，改革开放以来，我国农业生产关系和生产力水平发生了剧烈变化。使用上述的 SFA 方法不能很好反映政策变革、技术变化和投入变动对农业生产力的影响。

龚斌磊 2018 年在《 Journal of Development Economics 》上发表的《Agricultural reforms and production in China: Changes in provincial production function and productivity in 1978–2015》提出了一种新的效率分析方法——变系数 SFA（Varying Coefficient SFA） 。变系数 SFA 通过半参数方法，估计随时间和个体变化的 TFP 和要素弹性。

一、变系数 SFA 的简单介绍

1. 适用范围与特点

弹性生产函数。 传统的 SFA 假设生产函数是固定的，无法捕捉随时间变化的投入产出关系。变系数 SFA 允许生产函数的系数随时间变化，能够更好地反映农业生产中随着制度改革、技术进步等因素的变化。

控制外生冲击影响。 与传统 SFA 一样，由于使用 SFA 估计，变系数 SFA 仍然可以通过引入残差项控制外生冲击的影响，这种估计方程设定更加适合运用在农业的效率分析中。

半参数估计方法。 变系数 SFA 为两步估计，在估计变系数时，需要引入额外的变量来进行半参数估计。在这篇文章中，龚斌磊通过引入各地区农业、林业、牧业（渔业忽略）的增加值占比来估计变系数。

2. 分析内容

与传统 SFA 相同，变系数 SFA 可以估计要素产出弹性和技术效率，并进一步计算 TFP 和 TFP 增长率。

简单介绍一下变系数 SFA 的估计方程。Hastie and Tibshirani (1993) 首先引入了变系数模型（VCM），其中变系数是某些“阈值”变量 $\theta$ 的非参数函数。

其中，$\theta_1, \dots, \theta_p$ 通过未指定函数形式的 $h_1(\cdot), \dots, h_p(\cdot)$，改变 $x_1, \dots , x_p$ 系数。

其中：

  • $y_{it}$ 表示产出；
  • $x_{it}^k$ 表示第 $k$ 种投入要素；
  • $h_k(\theta_{it})$ 表示非参数函数，用于估计第 $k$ 种投入要素的变弹性 $\beta_{it}^k$。

二、变系数 SFA 的简明教程

我们重点关注怎么使用变系数 SFA。

1. 基础处理

该部分：（1）导入相关 R​ 包；（2）简单处理数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
library('gss')
library('Matrix')
library('np')
library('frontier')
library('plm')
library('mgcv')
library('gamlss')
library('assist')
library('sca')
library('AER')
library(bootBCa)

setwd("/Users/fgz/Desktop/变系数SFA/Code and Data (JDE)")  # 设置工作路径
t1=1978  # 开始时间
t2=2015  # 结束时间
data=read.csv("data.csv",header=T,fileEncoding = "GBK")  # 读取CSV数据
# 生成投入和产出的对数值
data$Y=log(data$GVAO/data$PPI/10)
data$Lr=log(data$labor/100)
data$Ld=log(data$land/1000)
data$Fr=log(data$fertiliser/100)
data$Kp=log(data$machinery/100)
# 生成新变量
data$total=data$farming+data$forestry+data$animalhusbandry+data$fisheries  # 农林牧渔总产值
data$ratio1=data$farming/data$total*100   # 农业占比
data$ratio2=data$forestry/data$total*100  # 林业占比
data$ratio3=data$animalhusbandry/data$total*100  #牧业占比
data$ratio4=data$fisheries/data$total*100  # 渔业占比

# 选定指定时间范围的数据，数据索引按照「行,列」进行，「,」表示数据检索只按照行检索，选择所有列
data=data[which(data$Year>=t1&data$Year<=t2),]
data$t=data$Year  # 生成一个新变量t, 取值等于year
second=data[,c(2:4,14:16)] # 选择指定变量
data=data[,c(2:4,18:29)]   # 选择指定变量
pdata=plm.data(data, c("SID","Year"))  # 创建面板数据, 指定ID和年份变量

2. 传统 SFA

控制时间固定效应的传统 SFA 代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
single <- sfa(Y~Lr+Ld+Fr+Kp+Year, timeEffect=TRUE, truncNorm=TRUE, data=pdata)
summary(single)
single_lr=coefficients(single)[2]
single_ld=coefficients(single)[3]
single_fr=coefficients(single)[4]
single_kp=coefficients(single)[5]
efficiencies(single)[,c(1,38)]
yy=t(efficiencies(single))
dim(yy) <- c(length(yy),1)
pdata$singleeff=yy
res=t(single$resid)
dim(res)=c(1178,1)
pdata$singleres=res

newdat=pdata
newdat$SID=11

newdat1=data
newdat1$Year=newdat1$Year+1
newdat1$t=newdat1$t+1
newdat1=plm.data(newdat1, c("SID","Year"))

3. 变系数 SFA

变系数 SFA 的第一步是要进行半参数估计，得到 $h_k(\theta_{it})$，这一步使用 GAMLSS​ 模型进行估计。第二步是基于 GAMLSS​ 模型的结果，计算变系数、残差、效率值等。由于 GAMLSS ​模型无法估计置信区间，需要使用 Bootstrap ​法进行估计。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# 定义模型：使用 Generalized Additive Models for Location, Scale and Shape (GAMLSS) 拟合模型
# 参考: https://psy-data.github.io/2015/theory_page6.html
# pvc: 变系数模型
model <- gamlss(Y~pvc(ratio1+ratio2+ratio3+t, by=Lr)
                +pvc(ratio1+ratio2+ratio3+t, by=Ld)
                +pvc(ratio1+ratio2+ratio3+t, by=Fr)
                +pvc(ratio1+ratio2+ratio3+t, by=Kp)
                +cs(ratio1+ratio2+ratio3+t)+Year+SID, data=pdata)

# 提取要素的边际贡献: 提取变系数估计结果, 除以中心化后的数值
pdata$lr=model$mu.s[,1]/(pdata$Lr-mean(pdata$Lr))
pdata$ld=model$mu.s[,2]/(pdata$Ld-mean(pdata$Ld))
pdata$fr=model$mu.s[,3]/(pdata$Fr-mean(pdata$Fr))
pdata$kp=model$mu.s[,4]/(pdata$Kp-mean(pdata$Kp))

# 计算估计残差, SFA 估计
# 两步估计有两个残差来源: gamlss、SFA
pdata$res1=fitted(model)-predict(model,newdata=newdat)  # gamlss残差
pdata$res1=pdata$res1+model$resid   # SFA残差
# SFA 估计
fit=sfa(res1 ~ 1, timeEffect=T, truncNorm=TRUE, data = pdata)
summary(fit)
# 提取 SFA 结果
# 提取效率值
yy=t(efficiencies(fit))
dim(yy) <- c(length(yy),1)
pdata$varyingeff=yy
# 提取残差
res=t(fit$resid)
dim(res)=c(1178,1)
pdata$varyingres=res

# Bootstrap 估计置信区间
# 创建两个矩阵，存储置信区间上限和下限
bootstrap_upper=matrix(0,4,(t2-t1+1))
bootstrap_lower=matrix(0,4,(t2-t1+1))
# 逐年估计 5% 置信区间, [5]-上限,[4]-下限
# delta = NA, 固定 bootstrap 次数, 而非自适应
for (i in t1:t2)
{
  # lr 置信区间上下限
  bootstrap_upper[1,(i-t1+1)]=BCa(pdata$lr[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[5]
  bootstrap_lower[1,(i-t1+1)]=BCa(pdata$lr[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[4]
  # ld 置信区间上下限
  bootstrap_upper[2,(i-t1+1)]=BCa(pdata$ld[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[5]
  bootstrap_lower[2,(i-t1+1)]=BCa(pdata$ld[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[4]
  # fr 置信区间上下限
  bootstrap_upper[3,(i-t1+1)]=BCa(pdata$fr[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[5]
  bootstrap_lower[3,(i-t1+1)]=BCa(pdata$fr[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[4]
  # kp 置信区间上下限
  bootstrap_upper[4,(i-t1+1)]=BCa(pdata$kp[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[5]
  bootstrap_lower[4,(i-t1+1)]=BCa(pdata$kp[which(pdata$Year==i)],NA,mean,alpha=c(0.025,0.975),M=10000)[4]
}

三、结果分析

上述代码已经完成变系数 SFA 的估计，但是由于模型为每个样本、每个要素都计算了一个系数，无法通过传统 SFA 的方式展示结果，需要通过图表的方式进行展示。下面是对上面结果的分析。

1. 年度要素弹性比较

该部分绘制了劳动、土地、化肥和机械四种要素的产出弹性年度变化图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
# Figure 2 Average Elasticities of the Four Inputs across Time
# 设置图形参数，创建 2x2 的多图布局
par(mfrow=c(2,2))
# 绘制劳动弹性图
plot(c(t1:t2),tapply(pdata$lr,pdata$Year,mean),type="l",xlab='Year',ylab='Labor Elasticity',lwd=3) # type="l" 绘制线条
# 绘制劳动弹性图: bootstrap 置信区间上限和下限
lines(c(t1:t2),bootstrap_upper[1,],type="l",lty=4,lwd=2) # 绘制虚线(点划线-lty=4)
lines(c(t1:t2),bootstrap_lower[1,],type="l",lty=4,lwd=2)
# 添加垂直参考线，用于标记年份
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)  # v-添加垂直线;h-添加水平线

# 绘制土地弹性图
plot(c(t1:t2),tapply(pdata$ld,pdata$Year,mean),type="l",xlab='Year',ylab='Land Elasticity',lwd=3,ylim=c(0.21,0.3))
# 绘制土地弹性图: bootstrap 置信区间上限和下限
lines(c(t1:t2),bootstrap_upper[2,],type="l",lty=4,lwd=2)
lines(c(t1:t2),bootstrap_lower[2,],type="l",lty=4,lwd=2)
# 添加垂直参考线，用于标记年份
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)

# 绘制化肥弹性图
plot(c(t1:t2),tapply(pdata$fr,pdata$Year,mean),type="l",xlab='Year',ylab='Fertiliser Elasticity',lwd=3)
# 绘制化肥弹性图: bootstrap 置信区间上限和下限
lines(c(t1:t2),bootstrap_upper[3,],type="l",lty=4,lwd=2)
lines(c(t1:t2),bootstrap_lower[3,],type="l",lty=4,lwd=2)
# 添加垂直参考线，用于标记年份
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)

# 绘制机械弹性图
plot(c(t1:t2),tapply(pdata$kp,pdata$Year,mean),type="l",xlab='Year',ylab='Machinery Elasticity',lwd=3)
# 绘制机械弹性图: bootstrap 置信区间上限和下限
lines(c(t1:t2),bootstrap_upper[4,],type="l",lty=4,lwd=2)
lines(c(t1:t2),bootstrap_lower[4,],type="l",lty=4,lwd=2)
# 添加垂直参考线，用于标记年份
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)

结果展示：

Fig 2: 年度要素弹性比较

2. 省际要素弹性比较

在完成年度比较后，可以从个体的维度（省份）去考察区域间要素产出弹性差异。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Figure 3 Average Elasticities of the Four Inputs across Provinces
state=c("北京","天津","河北","山西","内蒙","辽宁","吉林","黑龙江","上海","江苏","浙江","安徽","福建","江西","山东","河南","湖北","湖南","广东","广西","海南","重庆","四川","贵州","云南","西藏","陕西","甘肃","青海","宁夏","新疆")
# Create the input vectors.
# 设置图形参数，创建一个图形设备，并设置边距
par(mfrow = c(1, 1),mar=c(3, 10, 1, 3))
# 定义颜色向量
colors <- gray(0:4/4)[1:4]
# 定义输入向量
inputs <- c("Labor","Land","Fertilizer","Machinery")
# 创建一个值的矩阵，用于存储平均弹性和省份名称
Values <- matrix(0,5,31)
# 计算每个省份的平均弹性
Values[1,]=tapply(pdata$lr,pdata$SID,mean)
Values[2,]=tapply(pdata$ld,pdata$SID,mean)
Values[3,]=tapply(pdata$fr,pdata$SID,mean)
Values[4,]=tapply(pdata$kp,pdata$SID,mean)
# 设置省份名称
Values[5,]=c("Beijing","Tianjin","Hebei","Shanxi","Inner Mongolia","Liaoning","Jilin","Heilongjiang","Shanghai","Jiangsu","Zhejiang","Anhui","Fujian","Jiangxi","Shandong","Henan","Hubei","Hunan","Guangdong","Guangxi","Hainan","chongqing","Sichuan","Guizhou","Yunnan","Tibet"," Shaanxi","Gansu","Qinghai","Ningxia","Xinjiang")
# 调整列的顺序以匹配省份名称
Values=Values[,c(1,3,20,10,2,15,6,9,21,19,11,13,NA,4,5,7,8,16,18,12,14,17,NA,22:31)]
# 绘制条形图，水平放置
barplot(Values[1:4,],beside=F,horiz=T,las=1.5,names.arg=Values[5,],xlab="",ylab="",col=colors,xlim=c(0,1))
legend("bottomright", inputs, cex=1, fill=colors)  # 添加图例
dev.off()  # # 关闭图形设备

图片展示：

Fig 3: 省际要素弹性比较

3. 年度技术变化分析

上面两个部分都是针对要素产出弹性的分析，下面是对技术变化的年度分析。这里需要根据上面的结果，进一步计算技术变化指标。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# Technology Growth
par(mfrow = c(1, 1))
pdata$varyinggrowth[which(pdata$Year!=t1)]=exp(predict(model,newdata=newdat1[which(newdat1$Year!=(t2+1)),]))/exp(predict(model,newdata=pdata[which(pdata$Year!=t2),]))-1   #the growth rate of the next year
pdata$varyinggrowthtl[which(pdata$Year!=t1)]=exp(predict(model2,newdata=newdat1[which(newdat1$Year!=(t2+1)),]))/exp(predict(model2,newdata=pdata[which(pdata$Year!=t2),]))-1   #the growth rate of the next year
tg=c(0,coefficients(single)[6:(t2-t1+5)])
pdata$singlegrowth=NA
for (i in (t1+1):t2)
{
  pdata$singlegrowth[which(pdata$Year==i)]=tg[i-t1+1]-tg[i-t1]
}
growth=matrix(0,t2-t1+3,4)
growth[1:(t2-t1+1),1]=c(t1:t2)
growth[1:(t2-t1+1),2]=tapply(pdata$singlegrowth,pdata$Year,mean)
growth[1:(t2-t1+1),3]=tapply(pdata$varyinggrowth,pdata$Year,mean)

for (i in 2:38)
{
  growth[i,4]=sum(growth[2:i,3])
}
growth[(t2-t1)+2,2:3]=colSums(growth[1:(t2-t1+1),2:3],na.rm=T)
growth[(t2-t1)+3,2:3]=growth[(t2-t1+2),2:3]/(t2-t1)
growth  


# Figure 4 Technical Changes in the Six Reform Periods
plot(growth[1:38,1],exp(growth[1:38,4]),type="s",xaxt="n",yaxt="n",ylim=c(1,5))
axis(2, at=pretty(exp(growth[1:38,4])), lab=percent(pretty(exp(growth[1:38,4]))))
axis(1,at=seq(1980,t2,5))
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)

图片展示：

Fig 4: 年度技术变化分析

4. 省际技术效率比较

再来考察一下技术效率。

1
2
3
4
5
6
7
8
9
10
11
# Efficency
pdata$varyingeff
tapply(pdata$varyingeff,pdata$SID,mean)

# Figure 5 Ranking of the Average Technical Efficiency across Provinces 
par(mfrow = c(1, 1),mar=c(3, 10, 1, 3))
Values <- matrix(0,2,31)
Values[1,]=tapply(pdata$varyingeff,pdata$SID,mean)
Values[2,]=c("Beijing","Tianjin","Hebei","Shanxi","Inner Mongolia","Liaoning","Jilin","Heilongjiang","Shanghai","Jiangsu","Zhejiang","Anhui","Fujian","Jiangxi","Shandong","Henan","Hubei","Hunan","Guangdong","Guangxi","Hainan","chongqing","Sichuan","Guizhou","Yunnan","Tibet"," Shaanxi","Gansu","Qinghai","Ningxia","Xinjiang")
Values=Values[,order(Values[1,])]
barplot(as.numeric(Values[1,]),horiz=T,beside=F,las=1.5,names.arg=Values[2,],xlab="",ylab="",xlim=c(0,1))

图片展示：

Fig 5: 省际技术效率比较

5. TFP 增长率变化

再来看看 TFP 增长率的年度变化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# TFP
pdata$varyingtfp=exp(fitted(model)-pdata$lr*pdata$Lr-pdata$ld*pdata$Ld-pdata$fr*pdata$Fr-pdata$kp*pdata$Kp)
for (i in (t1+1):t2)
{
  pdata$varyingtfp[which(pdata$Year==i)]=pdata$varyingtfp[which(pdata$Year==(i-1))]/pdata$varyingeff[which(pdata$Year==(i-1))]*(pdata$varyinggrowth[which(pdata$Year==i)]+1)*pdata$varyingeff[which(pdata$Year==i)]
}
pdata$varyingtfpgrowth=NA
for (i in (t1+1):t2)
{
  pdata$varyingtfpgrowth[which(pdata$Year==i)]=pdata$varyingtfp[which(pdata$Year==i)]/pdata$varyingtfp[which(pdata$Year==(i-1))]-1
}
tapply(pdata$varyingtfpgrowth,pdata$Year,mean)

# Figure 6 Growth Rate of Total Factor Productivity across Time   
plot(c(t1:t2),tapply(pdata$varyingtfpgrowth,pdata$Year,mean), xaxt="n",yaxt="n",type = "h", lwd = 15, lend = 1)
abline(v=c(1984.5,1989.5,1993.5,1998.5,2003.5),lwd=1)
axis(2, at=pretty(tapply(pdata$varyingtfpgrowth,pdata$Year,mean)), lab=percent(pretty(tapply(pdata$varyingtfpgrowth,pdata$Year,mean))))
axis(1,at=seq(1980,t2,5))                   #figure 6

图片展示：

Fig 6：TFP增长率年度变化

四、写在后面

在这篇文章的后半部分，还基于上面的结果进行了实证分析。这里不在赘述，有兴趣可以查看原文。

下面附上文章 PDF、完整代码和数据，可用于代码复现。官方的代码有一些小bug，且没有代码备注，检验使用「Code and Data (JDE)」文件夹内的代码复现。

下载链接：百度网盘链接

‍