Stata：使用非线性最小二乘法估计服务需求强度函数模型

984 字 约 4 分钟 数据资料 Stata计量与统计方法数据定制与服务编程技巧

最近有小伙伴问了这样的一个问题：

她还给出了文献来源与示例数据。附件中的 pdf 文件就是出处：

A global analysis of residential heating and cooling service demand and cost-effective energy consumption under different climate change scenarios up to 2050

这个方程是该文章中提出的一个服务需求强度模型：

显然这是一个 SDI 关于 GDP 和 DD 的非线性方程。

在之前的系列课程「Stata 编程导论」的第 20 课时，我们讲解了使用 Stata 编写 nl 程序进行非线性最小二乘法的估计。这里我们就可以用上了！

×

×

首先，定义一个 nl 子程序：

cap prog drop nlsdi
prog def nlsdi
    version 14.0
    syntax varlist(numeric min = 3 max = 3) if, at(name)
    args y gpd hdd
    tempname C alpha beta gamma delta
    tempvar fterm gterm
    scalar `C' = `at'[1, 1]
    scalar `alpha' = `at'[1, 2]
    scalar `beta' = `at'[1, 3]
    scalar `gamma' = `at'[1, 4]
    scalar `delta' = `at'[1, 5]
    gen double `fterm' = 1 / (1 + exp(`alpha' * (`gdp' - `beta')))
    gen double `gterm' = (`hdd'^(`gamma') + `delta')
    *- 这里我就直接假设 M = 1 了
    replace `y' = `C' * `fterm' * `gterm' * 1
end

注意该程序的名称必须要以 nl 开头。这段代码基本是 nl 程序的模板，args y gpd hdd 里面的 y gpd hdd 就是要输入的三个变量，y 就是 SDI；C、alpha、beta、gamma 和 delta 是要估计的参数，这里为了思路更清晰，分别先计算 fterm（也就是 f(GDP) 函数）和 gterm（也就是 g(DD) 函数），然后在计算 y。

然后我们就可以读取数据使用这个子程序进行估计了：

import delimited using 数据.csv, clear
list in 1/10

*>     +-------------------------------------+
*>     | v1          y        gdp        hdd |
*>     |-------------------------------------|
*>  1. |  1   1954.333      20569   2257.789 |
*>  2. |  2   1115.376   18481.63    1995.11 |
*>  3. |  3   1340.796   11247.31   1639.053 |
*>  4. |  4   1266.477   10389.04   1853.882 |
*>  5. |  5   911.0911   13935.08    2152.04 |
*>     |-------------------------------------|
*>  6. |  6   1118.612   11094.61   1934.891 |
*>  7. |  7   1111.224   8718.498   1721.792 |
*>  8. |  8   862.0304   13625.73   2096.521 |
*>  9. |  9   836.1245      13491   2748.542 |
*> 10. | 10   1500.507   10514.06   2602.332 |
*>     +-------------------------------------+

*- 这里需要多次尝试更换初始值
nl sdi @ y gdp hdd, parameters(C alpha beta gamma delta) ///
    initial(C 10000 alpha -1 beta 20 gamma 0.5 delta 5) nolog

*>       Source |      SS            df       MS
*> -------------+----------------------------------    Number of obs =         46
*>        Model |  9.313e-10          0           .    R-squared     =     0.0000
*>     Residual |  5204207.8         45  115649.062    Adj R-squared =     0.0000
*> -------------+----------------------------------    Root MSE      =   340.0721
*>        Total |  5204207.8         45  115649.062    Res. dev.     =   665.8138

*> ------------------------------------------------------------------------------
*>            y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
*> -------------+----------------------------------------------------------------
*>           /C |      10000          .        .       .            .           .
*>       /alpha |  -.2251947   .0023205   -97.05   0.000    -.2298684   -.2205211
*>        /beta |         20          .        .       .            .           .
*>       /gamma |  -61.96606          .        .       .            .           .
*>       /delta |   9.980762          .        .       .            .           .
*> ------------------------------------------------------------------------------
*> Note: Parameter delta is used as a constant term during estimation.

*- 由于该模型的参数很多，提供的样本数据观测值又很少，所以该模型存在多个解
*- 使用不同的初始值可以得到不同的结果
nl sdi @ y gdp hdd, parameters(C alpha beta gamma delta) ///
    initial(C 20000 alpha -1 beta 20 gamma 0.5 delta 5) nolog

*>       Source |      SS            df       MS
*> -------------+----------------------------------    Number of obs =         46
*>        Model |  9.313e-10          0           .    R-squared     =     0.0000
*>     Residual |  5204207.8         45  115649.062    Adj R-squared =     0.0000
*> -------------+----------------------------------    Root MSE      =   340.0721
*>        Total |  5204207.8         45  115649.062    Res. dev.     =   665.8138
*>
*> ------------------------------------------------------------------------------
*>            y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
*> -------------+----------------------------------------------------------------
*>           /C |      20000          .        .       .            .           .
*>       /alpha |   -.260128   .0023077  -112.72   0.000    -.2647759   -.2554801
*>        /beta |         20          .        .       .            .           .
*>       /gamma |  -61.96658          .        .       .            .           .
*>       /delta |   9.980762          .        .       .            .           .
*> ------------------------------------------------------------------------------
*> Note: Parameter delta is used as a constant term during estimation.

nl sdi @ y gdp hdd, parameters(C alpha beta gamma delta) ///
    initial(C 20000 alpha -1 beta 30 gamma 0.5 delta 5) nolog

*>       Source |      SS            df       MS
*> -------------+----------------------------------    Number of obs =         46
*>        Model |  9.313e-10          0           .    R-squared     =     0.0000
*>     Residual |  5204207.8         45  115649.062    Adj R-squared =     0.0000
*> -------------+----------------------------------    Root MSE      =   340.0721
*>        Total |  5204207.8         45  115649.062    Res. dev.     =   665.8138
*>
*> ------------------------------------------------------------------------------
*>            y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
*> -------------+----------------------------------------------------------------
*>           /C |      20000          .        .       .            .           .
*>       /alpha |   -.178732   .0015372  -116.27   0.000    -.1818281   -.1756359
*>        /beta |         30          .        .       .            .           .
*>       /gamma |  -9114.889          .        .       .            .           .
*>       /delta |   11.69607          .        .       .            .           .
*> ------------------------------------------------------------------------------
*> Note: Parameter delta is used as a constant term during estimation.

其中 initial() 用来设定迭代的初始值，由于该模型的参数很多，提供的样本数据观测值又很少，所以该模型存在多个解，所以不同的初始值可能会拟合出不同的结果。实际估计中最好多做尝试。

点击这里跳转到 RStata 短书平台获取附件：Stata：使用非线性最小二乘法估计服务需求强度函数模型

相关文章

评论