Chapter2非线性最小二乘法与数值最优化

1、第1章非线性最小二乘法与数值最优化变量之间的关系更多地表现为非线性特征。线性模型作为基础模型是非线性的近似，即任何非线性模型都可以通过线性模型来近似表达。比如，模型y=B0+B1ex+u通过泰勒级数展开表述为y：oex|x田(x-xo)uo-ie"xo-,exox'u=：o*，：i*xu模型y=禺+Rx2+u的线性近似表达式为y:-o-i(2x)|xo(x-xo)u二o-2-1xo221xoxuo*，u例1.1利用MonteCarlo模拟的方法观察线性模型对非线性模型的近似。设DGP为：y=1O+O.2*exp(x)+u,x为1,3区间的均匀分布。利用线性模型与指数模型分别回

2、归模型，并计算x对y的(平均)边际影响与(平均)弹性。(数据文件：nonlin)但线性模型对非线性模型的近似程度取决于高阶部分是否充分小。即使在样本内线性模型能够较好地拟合数据，也不能准确地体现变量的结构关系。非线性模型中，x对y的边际影响(或弹性)是变化的；而线性模型中，x对y的边际影响(或弹性)是常数。很多情况下，线性模型与非线性模型对边际影响或弹性的估计存在非常大的差异。另外，利用线性模型拟合非线性数据存在潜在的危险，即区间外预测会存在越来越大的误差。因此，正确设定模型的形式是进行准确推断和预测的重要环节。对于一般的回归模型，如以下形式的模型，y=f(X,(3)+u1.1OLS一般不能得

3、到其解析解。比如，运用OLS方法估计模型(1.1),令S(P)表示残差平方和，即nn1.2S(四=2u：=£yi-f(Xi;B)2iWi1最小化S(P),即根据一阶条件可以得到；：S(B)郊9n开(Xi;=2JyT(Xi;0=。以模型y+Px'+u为例，其一阶条件为-S(3)1，“一2o-le2xi=0空份=一：e¥一昆一即函#=01-12i1乎y-0ie“2Xe2x=022i1上述方程组没有解析解，需要一般的最优化方法。很多数值最优化算法都可以完成这一类任务，这些方法的总体思路是一样的。即，从初始值出发，按照一定的方向搜寻更好的估计量，并反复迭代直至收敛。各种不同

4、的最优化算法的差异主要体现在三个方面：搜寻的方向、估计量变化的幅度和迭代停止法则。本章主要介绍非线性最小二乘法、梯度最优化方法及其其在Stata中的实现。1.1非线性最小二乘法非线性最小二乘法的思路是，通过泰勒级数将均值函数展开为线性模型。即，只包括一阶展开式，高阶展开式都归入误差项。然后再进行OLS回归，将得到的估计量作为新的展开点，再对线性部分进行估计。如此往复，直至收敛。设模型中共存在(k+1)个参数P=(P。，良，Pk)。首先为参数选择一组初始值甬=(凡,0邛1,0，Pk,o)。其中下标零表示初始值。然后将f(x,P)按泰勒级数在P。点展开。f(x,B)=f(X，自)+g(0)'

5、;(BM)+R1.3其中g0表示一阶导数在用=(?0,0,Pi,0，1,0)时的取值，R为高阶部分。上式中只保留P的线性部分，将高阶部分归入误差项，可得y=f(x,份+u=f(X，自)知丫”的)-Ru1.4=g(0)'B-f(x，隹)-g(0)M9其中，随机扰动项U1包含u和泰勒级数展开式中高阶成分。得到新的回归模型y-f(x,饱)'的)=g(0)pu11.5新的目标函数为Q*=?/«。模型的OLS估计量为，101)=g(0)g(0)'g(0)y_f(X，90),g(0),自0)工6=口0)'g(0)g(0)'g(0)y-f(X,00)因此，N

6、LS的迭代估计式为：1口j1)=瓦),g(j)g(j)1%)y-f(X，打)1.7上述方法也被称作高斯牛顿方法，即利用迭代方法估计非线性最小二乘模型。在(1.7)式中，第二项恰好是模型y-f(X,)=g(j)丫+v1.8丫的OLS估计量。把1.8式称作高斯-牛顿回归。P的迭代公式又可以写作口j书=瓦)+丫1.9这种迭代估计方法必须设定初始值和停止法则。初始值的选择对于迅速找到最优解非常重要。如果目标函数不是严格的凹函数或凸函数，或者存在多个局部最优值，可以设定多个初始值，观察最优解。如果不同的初始值得到相同的最优解，则结论是比较稳健的。停止法则用以设定满足一定的标准后终止迭代过程，否则迭代过程

7、会无限继续下去。可用的停止法则包括：目标函数Q(ftj+i)-Q(氏)没有明显的变化，或者g(j)的每个元素都非常小，或者氏+i)-氏)没有明显变化等。迭代法则也可以同时设定最高迭代次数n,如果经过n次迭代仍然没有能够达到收敛，则停止迭代。例1.2利用NLS方法估计非线性消费函数。(数据文件：usmacro)cs-:-incu利用NLS方法估计模型。.nl(realcons=alpha+beta*realgdpAgamma)参数初始值可以通过两种方法进行设定。其一，直接在参数表达式中设定。比如，上述模型中设定alpha=0,beta=0.5,gamma=1。Stata命令可以表述为：.nl(r

8、ealcons=alpha=0+beta=0.5*realgdpAgamma=1)另外一种方法是通过命令选项initial进行设定。比如上述命令也可以表述为：.nl(realcons=alpha+beta*realgdpAgamma),initial(alpha0beta0.5gamma1)注意，Stata的nl命令用于估计不存在缺失值的区间。如果某些变量存在缺失值，则必须通过variables选项设定模型中的解释变量。比如，.nl(realcons=alpha+beta*xAgamma=1),variables(realgdp)在一般的非线性模型中，解释变量对被解释变量的结构影响(边际影响或

9、弹性)都不是常数，甚至影响的方向与参数估计量的符号也不同。比如，在本模型中，收入对消费的边际影响为P/xinc；"Stata提供了一个方便的计算边际影响的命令mfx,用于计算不同变量在不同点上的边际影响。但如果在命令栏中使用nl命令，必须设定variables选项，才能使用mfx命令。比如，计算平均边际消费倾向.nl(realcons=alpha+beta*xAgamma=1),variables(realgdp).mfx,at(mean)计算收入inc=3000美元时的边际消费倾向，.mfxat(realgdp=3000)例1.3利用NLS方法估计如下生产函数模型。文件包含产出、资

10、本、劳动力等数据。(数据文件：production)(1) 广义CD(GCD)生产函数lny叩-011nK2lnLu实际上，GCD模型为线性模型，仍然可以采用nl命令进行估计，与regress命令得到完全相同的结果。.nl(lnout=-theta*out+beta0+beta1*lnk+beta2*lnl)(2) 不变替代弹性(CES)生产函数lny=：0-1/ln、K一,(1-、儿一,u.nl(lnout=beta0-1/rho=1*ln(delta=0.5*capitalA(-rho)+(1-delta)*laborA(-rho)资本劳动力的要素替代弹性为1/(1+rho)。.nlcom

11、(1/(1+_brho)(3) 广义CES生产函数lnyyuB01/：ln、.K-d(1、.)L“u.nl(lnout=beta0-1/rho=1*ln(delta=0.5*capitalA(-rho)+(1-delta)*laboS(-rho)1.2最优化方法概述非线性最小二乘法是利用泰勒级数展开的一种迭代估计方法，它是一般数值最优化的特例。对于一般的非线性的模型的其它估计方法，如极大似然法或广义矩方法等，都涉及到数值最优化的问题。接下来，我们介绍数值最优化方法的一般原理及其在Stata中的实现。设参数为P=(ft,P1,a)，目标函数为Q(P)o估计参数口使目标函数最优化(极大化或极小化)

12、。1.2.1 格点搜寻法格点搜寻法即令参数在其取值范围内取不同的数值，从中选择使目标函数最优的数值作为参数估计量。比如，模型中只有一个参数P,已知P介于0,1。我们可以令P从0逐渐变化到1,每次递增0.01,即0,0.01,0.02,，1。对于每一个取值计算目标函数，进而得到使目标函数最优的油计量。格点搜寻法只适用于参数个数比较少，而且取值范围具有一定的信息的模型。如果参数个数较多，而且取值范围较大的话，那么格点搜寻会耗费很长时间。比如，5个参数，每个参数取10个数值，那么需要计算105次。1.2.2 迭代方法另外一类最优化方法是迭代方法,为例，迭代方法的基本公式为：其中一类最常见的方法即是梯

13、度方法。)(j1)=Rd+A(j)9(j)以极大化问题1.10其中，9(j)为一阶导数梯度矩阵(gradient),元素为勿/虎。根据泰勒级数的线性展开式Q(份=Q(3(0)+9(0)'(3-0),第j次估计量表示为P(j+1),则目标函数Q(0(j+1)在M的展开式为：Q(瓦书)=Q(0j)+g(j)'(0j+)M)+R1.11其中，R表示被忽略的高阶项。将迭代公式代入可得Q(0M)Q(比)=g'a(j)g+R1.12对于极大化问题，如果A为正定矩阵且R比较小的话，则Q(P(j+i)>Q(P(j)。对于极小化问题，如果A为负定矩阵且R比较小的话，则Q(P(j+i

14、)<Q(P(j)。1.2.3 其它方法梯度方法要求目标函数充分平滑，可以计算其梯度向量。如果目标函数不存在梯度向量，则需要利用其它最优化方法。比如，最小绝对值方法的目标函数是，nny=f(X,0)+u,S(B)=Z|uJ=Z|y-f(X1|1 1i4可以利用线性规划方法。其它的更复杂的模型需要用到模拟退火法(annealling)或者遗传算法(genetic)。1.3 梯度方法1.3.1 牛顿(拉弗森)方法初始值为3(0)=(Po,o，Pi,0，久,0)。利用泰勒级数展开得到目标函数的二阶近似表达式：*1Q(份=Q(的)+g°)'(3-o)+2(3-E)'H(o

15、)(3-M)其中，g(0)为一阶导数梯度矩阵(gradient),元素为刃/谭j；H为二阶导数海塞矩阵(Hessian),元素为Q/cPi。如果函数为二次(或一次)线性函数，则泰勒级数展开式为精确展开。根据一阶条件*._._._s(B)/郊=%)+H(0)("的)=0可得P1的第一次估计量1%)=%)-H(0)g(0)如果Q*为凹函数，即H为正定矩阵，则闻)是最小值。如果Q*可以作为Q的较好的近似的话，P(i)也是?(目标函数Q的估计量)的较好的近似。Newton方法即是利用公式1fij+)=&j)-H(j)g(j)进行迭代计算。根据泰勒级数的线性展开式Q(3=Q(B(0)+

16、g(0)'(BM),目标函数Q(P(j+i)在即的展开式为:Q(B(»)=Q(B(j)+g(/(许书一0j)+R1.13其中，R表示被忽略的高阶项。将迭代公式代入可得1_Q(ftji)-Q(j)-P(j)Hg(j)R1.14对于极小化问题，如果H为正定矩阵且R比较小的话，则Q(岛+1)<Q('对于极大化问题，如果H为负定矩阵且R比较小的话，则Q(B(j+1)>Q(氏)。实践中，梯度方法对应两种修正方法。其一，根据(1.5)式，如果A比较小，则目标函数的变化也比较小，每次迭代的进展速度会比较慢。反之，如果A比较大，则有可能导致高阶项R比较大，将其忽略会产生较

17、大的误差。因此，对梯度方法的一个修正是加上一个调整因子ao这种方法称之为修正的牛顿方法。1口j书=ftj)-a(j)H(j)qj)1.15建得下式最大化，Q?(a)=Q(M-aH(7)g(j)卜其二，如果出现H为奇异矩阵的情况，则上式无法进行计算。这时，需要对奇异矩阵进行调整。比如，令H加上某个矩阵H,得到D=(H+C)-1。利用下式进行迭代估计。1比+)=自j)-%)D(j)g(j)16这种方法称之为拟牛顿方法。其中，D(j)作为H(j)的近似，总是正定矩阵。如果D(j)=H(j),则得到了修正的Newton方法。如果D二H，W)=1,则得到了普通的Newton方法。多情况下，比如Q*不是严

18、格的凹函数，Newton方法可能不能正确地找到P的最优值。实际上，对比NR公式可以发现，高斯牛顿方法是牛顿拉弗森方法的一个特例。拟牛顿方法的基本步骤如下。第一步，计算g(j)和D(j),用于决定下一步的搜寻方向H(；g(j)O第二步，计算他)，常用计算方法是最小化如下函数?1Q(=Q(瓦)7H(j)g(j)第三步，根据停止法则，判断继续迭代还是停止迭代。1.3.2 得分法与NR方法相比，得分法将其中的海塞矩阵H替换为H的期望矩阵的估计量，即Hms1.3.3 BHHH法Berndt,Hall,Hall,andHausman(1974)方法将海塞矩阵H替换为HHBHHH?那c6pBHHH估计只需要

19、计算目标函数的一阶导数，而无需计算二阶导数，因此计算比较方便。1.3.4 陡峭爬山法陡峭爬山法直接令H=I。1.3.5 DFP、BFGS方法Davidson,Fletcher,andPowell是通过下面的迭代公式计算H矩阵,AA,AjYjYj'Aj'Aj1AjgjSj'YjYj'AjTj石二Aj-gj,二gj1-gjBoyden,Fletcher,Goldfarb,andShannon对其进行了进行了进一步修正。今今'AjTjTj'Aj'Aj*=Aj+”+gjjjjj-VAjTj为犷BjYj丫jAjYj.=(Bj/BjYj)-(AjjY

20、j'Aj%)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设定了三种类型的估值

21、器：v0、v1、v2o如果采用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中，并将

22、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",&qu

23、ot;v1","v2","v1debug","v2debug"。默认选项为d0。注：对于数学函数的最优化，evaluatortype包括"d0","d1","d2","d1debug","d2debug"。4 .最优化方法命令格式为：optimize_init_technique(S,technique)technique包括："nr"修正的Newton-Raphson"dfp"Davi

24、don-Fletcher-Powell"bfgs"Broyden-Fletcher-Goldfarb-Shanno"bhhh"Berndt-Hall-Hall-Hausman"nm"Nelder-Mead5 .奇异海塞矩阵情况的处理方法命令格式为：optimize_init_singularHmethod(S,singularHmethod)singularHmethod包括："m-marquardt"修正的Marquardt算法"hybrid"垂直爬升法和牛顿方法的混合6 .初始值的设定命令格

25、式为：optimize_init_params(S,realrowvectorinitialvalues)7 .停止法则的设定设置最高迭代次数：optimize_init_conv_maxiter(S,realscalarmax)设置收敛标准:对于nr、dfp、bfgs、bhhh方法，收敛标准为:(1)根据参数估计量的变化。即,optimize_init_conv_ptol(S,(2)根据目标函数的变化。即,optimize_init_conv_vtol(S,(3)根据一阶导数的变化。即,optimize_init_conv_nrtol(S,(4)对于极大化问题，要求大化问题。当满足条件(1)

26、、(3)、(4)默认的收敛标准为ptol=1e-68.设置跟踪程序命令格式为：optimize_init_tracelevel(S,其中，tracelevel包括：-Hmreldif(p,p_prior)<ptolrealscalarptol)reldif(v,v_prior)<vtol。realscalarvtol)g*invsym(-H)*g'<nrtol。realscalarvtol)为正半定矩阵。对于极小化问题，可以通过-f转为极时，或者满足(2)、(3)、(4)时，Stata将停止迭代。,vtol=1e-7,nrtol=1e-5stringscalartra

27、celevel)"none""value""tolerance""params""step""gradient""hessian"不进行记录目标函数值(默认选项)收敛值参数倍计值变化幅度梯度向量海塞矩阵1.4.3最优化并观察估计结果进行最优化的命令格式：p=optimize(S)Stata通过optimize_result_*观察估计结果，包括参数估计量及其协方差矩阵、目标函数值、梯度向量、得分向量、海塞矩阵、迭代次数、是否收敛等。参数估计量：realro

28、wvectoroptimize_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方法有意义)：rea

29、lmatrixoptimize_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,其

30、中S表示得分向量，即稳健协方差矩阵。迭代次数：realscalaroptimize_result_iterations(S)是否收敛：realscalaroptimize_result_converged(S)如果达到收敛标准，贝Uoptimize_result_converged()=1。观察初始设定：optimize_query(S)观察当前的收敛标准的设定值：realscalaroptimize_init_conv_ptol(S)realscalaroptimize_init_conv_vtol(S)realscalaroptimize_init_conv_nrtol(S)1.5最优化案

31、例例1.4最优化如下函数。miny-x2x-3一阶导数为g=-2x+1,二阶导数为H=-2O根据NR迭代公式比«=Wj)_H乙g(j),估计量的迭代计算公式为j书=x(j)-0.5x(-2x+1)|x(j)。设定初始值为x(0)=0,则x(i)0.5(-2x1)|他00.5=0.5x=)+0.5乂(一2x+1)(=0.5+0=0.5目标函数不再变化，停止迭代。事实上，由于目标函数为严格凸函数，初始点无论选择在哪里，都可以经过一次迭代达到收敛。Stata的最优化命令如下。(程序文件：optim_func1.do)。mata:mataclearvoidmyeval(todo,x,y,g,

32、H)(y=-xA2+x-3if(todo>=1)g=-2*x+1if(todo=2)H=-2)S=optimize_init()optimize_init_evaluator(S,&myeval()optimize_init_params(S,10)optimize_init_tracelevel(S,"step")x=optimize(S)end例1.5最优化如下函数。22Qtexp(x1x2x1x2-x1x23)梯度向量为fQ/%(2XiX2-1)Qg=J=£Q/%_(2x2+x+1)Q_海塞矩阵为cQcQ>18；：2q次1况!.心=,;：

33、2Q2Q(2xx2-1)2QQ(2x2X1)(2xx2-1)QQ(2x1x2-1)(2x2Xi1)Q2Q(2x2Xi1)2Qx2x2根据NR迭代公式，的加=8j)-Hj)g，估计量的迭代计算公式为设定初始值为Xi(0)=0,X2(0)=0。Stata的命令函数如下(程序文件:mata:optim_func2.do)。mataclearvoidmyeval(todo,x,y,g,H)(y=x1A2+x2A2+x1*x2卜x1+x2+3if(todo>=1)g1=(2*x1+x2卜1)*yg2=(2*x2+x1+1)*yif(todo=2)H1,1=2*y+(2*x1+x2-1)*g1H1,

34、2=y+(2*x2+x1+1)*g1H2,2=2*y+(2*x2+x1+1)*g2_makesymmetric(H)S=optimize_init()optimize_init_evaluator(S,&myeval()optimize_init_which(S,"min")optimize_init_params(S,(0,0)optimize_init_tracelevel(S,"value")x=optimize(S)end例1.6利用Stata的最优化方法估计例1.1中的消费函数模型。(数据文件:文件：usmacro.do)/*usmac

35、ro.dta*/mata:mataclearData=st_data(.,("realcons","realgdp")y=Data.,1e=J(rows(Data),1,1)X=e,Data.,2:cols(Data)voideval(todo,p,X,y,v,g,H)(beta0=p1beta1=p2beta2=p3mu=beta0*X.,1:+beta1*X.,2:Abeta2v=(y:-mu):A2)S=optimize_init()optimize_init_evaluator(S,&eval()optimize_init_which(S,&quo