Chapter2 非线性最小二乘法与数值最优化文档格式.docx
《Chapter2 非线性最小二乘法与数值最优化文档格式.docx》由会员分享，可在线阅读，更多相关《Chapter2 非线性最小二乘法与数值最优化文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
。
其中下标零表示初始值。
然后将f（X,）按泰勒级数在0点展开。
1.3
其中g0表示一阶导数在0=（0,0,1,0,…,k,0）时的取值，R为高阶部分。
上式中只保留的线性部分，将高阶部分归入误差项，可得
1.4
其中，随机扰动项u1包含u和泰勒级数展开式中高阶成分。
得到新的回归模型
1.5
新的目标函数为
模型的OLS估计量为，
1.6
因此，NLS的迭代估计式为：
1.7
上述方法也被称作高斯牛顿方法，即利用迭代方法估计非线性最小二乘模型。
在（1.7）式中，第二项恰好是模型
1.8
的OLS估计量。
把1.8式称作高斯-牛顿回归。
的迭代公式又可以写作
1.9
这种迭代估计方法必须设定初始值和停止法则。
初始值的选择对于迅速找到最优解非常重要。
如果目标函数不是严格的凹函数或凸函数，或者存在多个局部最优值，可以设定多个初始值，观察最优解。
如果不同的初始值得到相同的最优解，则结论是比较稳健的。
停止法则用以设定满足一定的标准后终止迭代过程，否则迭代过程会无限继续下去。
可用的停止法则包括：
目标函数Q（（j+1））-Q（（j））没有明显的变化，或者g（j）的每个元素都非常小，或者（j+1）-（j）没有明显变化等。
迭代法则也可以同时设定最高迭代次数n，如果经过n次迭代仍然没有能够达到收敛，则停止迭代。
例1.1利用NLS方法估计非线性消费函数。
usmacro）
利用NLS方法估计模型。
.nl（realcons={alpha}+{beta}*realgdp^{gamma}）
参数初始值可以通过两种方法进行设定。
其一，直接在参数表达式中设定。
比如，上述模型中设定alpha=0，beta=0.5，gamma=1。
Stata命令可以表述为：
.nl（realcons={alpha=0}+{beta=0.5}*realgdp^{gamma=1}）
另外一种方法是通过命令选项initial进行设定。
比如上述命令也可以表述为：
.nl（realcons={alpha}+{beta}*realgdp^{gamma}）,initial（alpha0beta0.5gamma1）
注意，Stata的nl命令用于估计不存在缺失值的区间。
如果某些变量存在缺失值，则必须通过variables选项设定模型中的解释变量。
比如，
.nl（realcons={alpha}+{beta}*x^{gamma=1}）,variables（realgdp）
在一般的非线性模型中，解释变量对被解释变量的结构影响（边际影响或弹性）都不是常数，甚至影响的方向与参数估计量的符号也不同。
比如，在本模型中，收入对消费的边际影响为inc-1。
Stata提供了一个方便的计算边际影响的命令mfx，用于计算不同变量在不同点上的边际影响。
但如果在命令栏中使用nl命令，必须设定variables选项，才能使用mfx命令。
比如，计算平均边际消费倾向
.mfx,at（mean）
计算收入inc=3000美元时的边际消费倾向，
.mfxat（realgdp=3000）
例1.2利用NLS方法估计如下生产函数模型。
文件包含产出、资本、劳动力等数据。
production）
（1）广义CD（GCD）生产函数
实际上，GCD模型为线性模型，仍然可以采用nl命令进行估计，与regress命令得到完全相同的结果。
.nl（lnout=-{theta}*out+{beta0}+{beta1}*lnk+{beta2}*lnl）
（2）不变替代弹性（CES）生产函数
.nl（lnout={beta0}-1/{rho=1}*ln（{delta=0.5}*capital^（-{rho}）+（1-{delta}）*labor^（-{rho}）））
资本劳动力的要素替代弹性为1/（1+rho）。
.nlcom（1/（1+_b[rho]））
（3）广义CES生产函数
1.2最优化方法概述
非线性最小二乘法是利用泰勒级数展开的一种迭代估计方法，它是一般数值最优化的特例。
对于一般的非线性的模型的其它估计方法，如极大似然法或广义矩方法等，都涉及到数值最优化的问题。
接下来，我们介绍数值最优化方法的一般原理及其在Stata中的实现。
设参数为=（0,1,…,k），目标函数为Q（）。
估计参数使目标函数最优化（极大化或极小化）。
1.2.1格点搜寻法
格点搜寻法即令参数在其取值范围内取不同的数值，从中选择使目标函数最优的数值作为参数估计量。
比如，模型中只有一个参数，已知介于[0，1]。
我们可以令从0逐渐变化到1，每次递增0.01，即0，0.01，0.02，…,1。
对于每一个取值计算目标函数，进而得到使目标函数最优的估计量。
格点搜寻法只适用于参数个数比较少，而且取值范围具有一定的信息的模型。
如果参数个数较多，而且取值范围较大的话，那么格点搜寻会耗费很长时间。
比如，5个参数，每个参数取10个数值，那么需要计算105次。
1.2.2迭代方法
另外一类最优化方法是迭代方法，其中一类最常见的方法即是梯度方法。
以极大化问题为例，迭代方法的基本公式为：
1.10
其中，g（j）为一阶导数梯度矩阵（gradient），元素为
根据泰勒级数的线性展开式
，第j次估计量表示为（j+1），则目标函数Q（（j+1））在（j）的展开式为：
1.11
其中，R表示被忽略的高阶项。
将迭代公式代入可得
1.12
对于极大化问题，如果A为正定矩阵且R比较小的话，则Q（（j+1））>
Q（（j））。
对于极小化问题，如果A为负定矩阵且R比较小的话，则Q（（j+1））<
1.2.3其它方法
梯度方法要求目标函数充分平滑，可以计算其梯度向量。
如果目标函数不存在梯度向量，则需要利用其它最优化方法。
比如，最小绝对值方法的目标函数是，
可以利用线性规划方法。
其它的更复杂的模型需要用到模拟退火法（annealling）或者遗传算法（genetic）。
1.3梯度方法
1.3.1牛顿（拉弗森）方法
初始值为
利用泰勒级数展开得到目标函数的二阶近似表达式：
其中，g（0）为一阶导数梯度矩阵（gradient），元素为
；
H（0）为二阶导数海塞矩阵（Hessian），元素为
如果函数为二次（或一次）线性函数，则泰勒级数展开式为精确展开。
根据一阶条件
可得的第一次估计量
如果Q*为凹函数，即H为正定矩阵，则
（1）是最小值。
如果Q*可以作为Q的较好的近似的话，
（1）也是
（目标函数Q的估计量）的较好的近似。
Newton方法即是利用公式
进行迭代计算。
，目标函数Q（（j+1））在（j）的展开式为：
1.13
1.14
对于极小化问题，如果H为正定矩阵且R比较小的话，则Q（（j+1））<
对于极大化问题，如果H为负定矩阵且R比较小的话，则Q（（j+1））>
实践中，梯度方法对应两种修正方法。
其一，根据（1.5）式，如果A比较小，则目标函数的变化也比较小，每次迭代的进展速度会比较慢。
反之，如果A比较大，则有可能导致高阶项R比较大，将其忽略会产生较大的误差。
因此，对梯度方法的一个修正是加上一个调整因子。
这种方法称之为修正的牛顿方法。
1.15
使得下式最大化，
其二，如果出现H为奇异矩阵的情况，则上式无法进行计算。
这时，需要对奇异矩阵进行调整。
比如，令H加上某个矩阵H，得到D=（H+C）-1。
利用下式进行迭代估计。
1.16
这种方法称之为拟牛顿方法。
其中，D（j）作为H（j）的近似，总是正定矩阵。
如果D（j）=H（j），则得到了修正的Newton方法。
如果D（j）=H（j），（j）=1，则得到了普通的Newton方法。
多情况下，比如Q*不是严格的凹函数，Newton方法可能不能正确地找到的最优值。
实际上，对比NR公式可以发现，高斯牛顿方法是牛顿拉弗森方法的一个特例。
拟牛顿方法的基本步骤如下。
第一步，计算g（j）和D（j），用于决定下一步的搜寻方向
第二步，计算（j），常用计算方法是最小化如下函数
第三步，根据停止法则，判断继续迭代还是停止迭代。
1.3.2得分法
与NR方法相比，得分法将其中的海塞矩阵H替换为H的期望矩阵的估计量，即
1.3.3BHHH法
Berndt,Hall,Hall,andHausman（1974）方法将海塞矩阵H替换为
BHHH估计只需要计算目标函数的一阶导数，而无需计算二阶导数，因此计算比较方便。
1.3.4陡峭爬山法
陡峭爬山法直接令H=I。
1.3.5DFP、BFGS方法
Davidson,Fletcher,andPowell是通过下面的迭代公式计算H矩阵，
Boyden,Fletcher,Goldfarb,andShannon对其进行了进行了进一步修正。
1.4最优化方法在Stata中的实现
Stata中通过Mata的最优化函数optimize对目标函数进行最优化。
基本步骤包括三步。
设定目标函数，设置最优过程，进行最优化并观察估计结果。
1.4.1设定目标函数
目标函数的一般格式为：
voidevaluator（todo,p,X1[,X2,...,Xk],v,g,H）
其中，todo为0、1、2。
0表示只设定目标函数（默认选项），1表示设定目标函数和梯度向量，2表示设定目标函数、梯度向量和海塞矩阵。
P为待估参数，v为目标函数。
g、H表示梯度向量和海塞矩阵。
X1[,X2,...,Xk]表示模型中的变量，最多可以设置9个。
Stata设定了三种类型的估值器：
v0、v1、v2。
如果采用v1，则需要用户设定梯度向量的表达式；
如果采用v2，则需要用户设定海塞矩阵的表达式。
模型中的变量必须通过下面的命令进行设置。
如果有两组变量X、y，则进行如下设置。
voidevaluator（todo,p,X,y,v,g,H）
optimize_init_argument（S,1,X）
optimize_init_argument（S,2,y）
1.4.2设置最优化过程
最优化过程的控制主要包括最优化类型（最大化、最小化）、估值器类型、最优化方法、初始值、迭代停止法则等。
1．初始化最优设置
命令格式为：
S=optimize_init（）
这一命令将最优化结果保存到S中，并将optimize_init_*（）的设置恢复到默认状态。
如果需要更改设置，通过后面的各种optimize_init_*（）进行单独设定。
optimize_init_evaluator（S,&
evaluator（））
这一命令将最优化设定赋值到目标函数
2．设置最优化类型
optimize_init_which（S,{"
max"
|"
min"
}）
默认选项为max。
3．设置估值器类型
optimize_init_evaluatortype（S,evaluatortype）
evaluatortype包括：
"
v0"
"
v1"
v2"
v1debug"
v2debug"
默认选项为d0。
注：
对于数学函数的最优化，evaluatortype包括"
d0"
d1"
d2"
d1debug"
d2debug"
4．最优化方法
optimize_init_technique（S[,technique]）
technique包括：
nr"
修正的Newton-Raphson
dfp"
Davidon-Fletcher-Powell
bfgs"
Broyden-Fletcher-Goldfarb-Shanno
bhhh"
Berndt-Hall-Hall-Hausman
nm"
Nelder-Mead
5．奇异海塞矩阵情况的处理方法
optimize_init_singularHmethod（S[,singularHmethod]）
singularHmethod包括：
m-marquardt"
修正的Marquardt算法
hybrid"
垂直爬升法和牛顿方法的混合
6．初始值的设定
optimize_init_params（S,realrowvectorinitialvalues）
7．停止法则的设定
设置最高迭代次数：
optimize_init_conv_maxiter（S,realscalarmax）
设置收敛标准：
对于nr、dfp、bfgs、bhhh方法，收敛标准为：
（1）根据参数估计量的变化。
即，mreldif（p,p_prior）<
ptol
optimize_init_conv_ptol（S,realscalarptol）
（2）根据目标函数的变化。
即，reldif（v,v_prior）<
vtol。
optimize_init_conv_vtol（S,realscalarvtol）
（3）根据一阶导数的变化。
即，g*invsym（-H）*g'
<
nrtol。
optimize_init_conv_nrtol（S,realscalarvtol）
（4）对于极大化问题，要求-H为正半定矩阵。
对于极小化问题，可以通过-f转为极大化问题。
当满足条件
（1）、（3）、（4）时，或者满足
（2）、（3）、（4）时，Stata将停止迭代。
默认的收敛标准为ptol=1e-6，vtol=1e-7，nrtol=1e-5。
8．设置跟踪程序
optimize_init_tracelevel（S,stringscalartracelevel）
其中，tracelevel包括：
------------------------------------------------------
"
none"
不进行记录
value"
目标函数值（默认选项）
tolerance"
收敛值
params"
参数估计值
step"
变化幅度
gradient"
梯度向量
hessian"
海塞矩阵
-----------------------------------------------------
1.4.3最优化并观察估计结果
进行最优化的命令格式：
p=optimize（S）
Stata通过optimize_result_*观察估计结果，包括参数估计量及其协方差矩阵、目标函数值、梯度向量、得分向量、海塞矩阵、迭代次数、是否收敛等。
参数估计量：
realrowvectoroptimize_result_params（S）
目标函数最优值：
realscalaroptimize_result_value（S）
目标函数初始值：
realscalarroptimize_result_value0（S）
梯度向量：
realrowvectoroptimize_result_gradient（S）
得分向量：
realmatrixoptimize_result_scores（S）
海塞矩阵：
realmatrixoptimize_result_Hessian（S）
协方差矩阵：
realmatrixoptimize_result_V（S）
协方差矩阵（仅对ML方法有意义）：
realmatrixoptimize_result_V_oim（S）
realmatrixoptimize_result_V_opg（S）
realmatrixoptimize_result_V_robust（S）
其中，oim表示Cov=invsym（-H）（最大化）（对于最小化问题，Cov=invsym（H）），即直接利用观测到的信息矩阵（observedinformatinmatrix）来计算。
Opg表示Cov=invsym（S'
S），其中S表示得分向量，即利用梯度向量的外积公式（outerproductofgradient）。
Robust表示Cov=H*invsym（S'
S）*H，其中S表示得分向量，即稳健协方差矩阵。
迭代次数：
realscalaroptimize_result_iterations（S）
是否收敛：
realscalaroptimize_result_converged（S）
如果达到收敛标准，则optimize_result_converged（）=1。
观察初始设定：
optimize_query（S）
观察当前的收敛标准的设定值：
realscalaroptimize_init_conv_ptol（S）
realscalaroptimize_init_conv_vtol（S）
realscalaroptimize_init_conv_nrtol（S）
1.5最优化案例
例1.3最优化如下函数。
一阶导数为g=-2x+1，二阶导数为H=-2。
根据NR迭代公式
，估计量的迭代计算公式为
设定初始值为x（0）=0，则
目标函数不再变化，停止迭代。
事实上，由于目标函数为严格凸函数，初始点无论选择在哪里，都可以经过一次迭代达到收敛。
Stata的最优化命令如下。
（程序文件：
optim_func1.do）。
mata:
mataclear
voidmyeval（todo,x,y,g,H）
{
y=-x^2+x-3
if（todo>
=1）{
g=-2*x+1
if（todo==2）{
H=-2
}
}
}
optimize_init_evaluator（S,&
myeval（））
optimize_init_params（S,10）
optimize_init_tracelevel（S,"
）
x=optimize（S）
end
例1.4最优化如下函数。
梯度向量为
海塞矩阵为
根据NR迭代公式，
设定初始值为x1（0）=0，x2（0）=0。
Stata的命令函数如下（程序文件：
optim_func2.do）。
y=x[1]^2+x[2]^2+x[1]*x[2]-x[1]+x[2]+3
g[1]=（2*x[1]+x[2]-1）*y
g[2]=（2*x[2]+x[1]+1）*y
H[1,1]=2*y+（2*x[1]+x[2]-1）*g[1]
H[1,2]=y+（2*x[2]+x[1]+1）*g[1]
H[2,2]=2*y+（2*x[2]+x[1]+1）*g[2]
_makesymmetric（H）
optimize_init_which（S,"
optimize_init_params（S,（0,0））
例1.5利用Stata的最优化方法估计例1.1中的消费函数模型。
usmacro；
程序文件：
usmacro.do）
/*usmacro.dta*/
mata:
mataclear
Data=st_data（.,（"
realcons"
realgdp"
））
y=Data[.,1]
e=J（rows（Data）,1,1）
X=e,Data[.,2:
:
cols（Data）]
voideval（todo,p,X,y,v,g,H）
{
beta0=p[1]
beta1=p[2]
beta2=p[3]
升级会员 