R语言实现随机前沿分析SFA、数据包络分析DEA、自由处置包分析FDH和BOOTSTRAP方法
原文链接：http://tecdat.cn/?p=23709
在经济学中，技术效率是指在既定的投入下产出可增加的能力或在既定的产出下投入可减少的能力。常用度量技术效率的方法是生产前沿分析方法。所谓生产前沿是指在一定的技术水平下，各种比例投入所对应的最大产出集合。而生产前沿通常用生产函数表示。前沿分析方法根据是否已知生产函数的具体的形式分为参数方法和非参数方法，前者以随机前沿分析（StochasticFrontierAnalysis，下文简称SFA）为代表，后者以数据包络分析（DataEnvelopeAnalysis，下文简称DEA）为代表。
本文在R软件中实现SFA、DEA与自由处置包分析FDH模型。
用R进行模拟
1.按照p=0.2的二项分布，随机生成一个大小为100的向量。
> x
plot(table(x), main = "frequency")
其他分布。泊松P（λ）（函数rpois）等。
2.数字变量 按照高斯分布N（μ=1，σ=1）随机生成一个大小为100的向量。
> x hist(x, main = "")
其他分布。Uniform U[a,b]（函数runif）等等。
模拟数据
输出y和输入x都是n=15的大小
真正的前沿是由函数定义的。
为了模拟数据:
1. 定义输入的矢量为x∼U[0,1]
2. 定义一个向量u∼N+(µ = 0.25, σ = 0.2)
3. 输入的矢量定义为
函数set.seed允许我们保持相同的模拟数据
> u = rtmvnorm(n = ns, mean = c(0.25), sigma = c(0.2),
+ lower = c(0))
> y = ybar/(1 + u)
数据的表示方法
模拟数据。
> plot(y ~ x
绘制真正的边界。
> lines(t.fr ~ x.seq)
"真实前沿 "效率测算
以产出为导向的测算。
输入导向的方法:
Shepard 方法:
> lambda = y/sqrt(x)
> theta = y^2/x
> delta = 1/theta
可重复研究
> matable <- xtable(tab1\[1:5, \], digits = 3, align = "l|ccc",
+ caption = ""真实前沿效率测算")
表:真实_前沿效率_
随机前沿
1.用函数lm调整一个线性模型，并保留回归线的系数β：y = α + βx
2. 找到能使（yi-yˆi）最大化的公司k，i=1，...，n。注意，公司k可以通过函数识别来手动找到和检测
3. 计算α 0，使回归线y = α 0 + βx穿过企业k，并代表随机前沿。
点击标题查阅相关内容
左右滑动查看更多
01
02
03
04
1. OLS 模型
> res.lm <- lm(y ~ x)
2. 使用函数识别
> identify(x, y)
3. 找到α0
> alpha2 <- y\[3\] - beta.lm\[2\] *
+ x\[3\]
> plot(y ~ x, type = "p",
+ 1))
随机前沿效率测算
让我们定义f1 : x → α 0 + βx
> f1 = function(x) alpha2 + beta.lm\[2\] * x
> f1.inv = function(x) (x - alpha2)/beta.lm\[2\]
面向产出的方法。
面向输入的方法:
Shepard 方法:
DEA - FDH 表示
手动检测位于两个边界上的公司
> plot(y ~ x)
> lines(x\[c(2, 9, 3, 4)\], legend = c("DEA",
+ "FDH"))
DEA-FDH效率前沿/衡量标准
让我们考虑5号公司
1. 如果该公司在输出方向上是有效的，它将位于前沿线的哪一部分？在输入方向上？
2. 利用这个估计前沿的位置，计算出效率的衡量标准
Naive Bootstrap
重复B次(用循环的方式)
1.用函数样本在15个观测值中取样
2.计算前沿的新估计值
3.计算新的效率方法
4. 储存结果，计算偏差, 方差, 置信区间
对真实数据的分析
数据集
在62个农场观察到一个输出变量和三个输入变量
> plot(y ~ x1,pch = 16, col = "blue")
> abline(lm(y ~ x1, data), col = "red")
散点图3-d
> scatterplot3d(x1,
+ x2, y)
数据结构
p个输入包含在一个p×n的矩阵中
> input <- t(cbind(spain$x1, spain$x2, spain$x3))
q个输出包含在一个q×n的矩阵中
> output <- t(matrix(spain$y))
效率的衡量标准
计算 DEA 效率估计值
计算 FDH 效率估计值
计算m阶效率估计值
计算非参数的条件和非条件的α-量化估计（默认情况下，α=0.95
> res <- rbind(dea, fdh, orderm,
+ res.hquan)
你可以使用函数order或sort来计算企业的排名，排名根据效率测算。
> plot(density(res.dea)
Bootstrap
函数boot实现了Simar和Wilson（1998）的bootstrap方法，用于估计Shepard（1970）输入和输出距离函数的置信区间。
> boot(input, output)