原文：https://bashtage.github.io/linearmodels/panel/index.html介绍面板数据是具有时间序列（T）和截面（N）两个维度组织的数据。在面板数据的大多数经典应用中，观测的数量N大，而时间段的数量T小（通常在2到5之间）。这些估计量的渐近理论多是在N固定而T会发散的前提下得到的的。大多数面板模型旨在估计下面方程的参数：$$ y_{i t}=x_{i t} \beta+\alpha_{i}+\epsilon_{i t} $$其中i为观测，t为时间，$\beta$包含了要估计的参数。$\alpha_i$是个体的异质性，无法观测。$\epsilon_{i t}$是特异性误差（idiosyncratic errors），与$\alpha_i$和$x_{it}$不相关。所有的模型都需要两个输入：dependent，即因变量$y_{it}$exog，即自变量$x_{it}$并使用不同的技术处理$\alpha_i$PanelOLS固定效应，即面板数据中随个体变化但不随时间变化的一类变量方法。这在数学上等价于为每个观测包含一个虚拟变量，尽管出于性能原因，实际上并不是这样实现的。BetweenOLS平均数估计。即对面板数据的每个个体求平均数，然后利用自变量和因变量的N组观测值估计参数。FirstDifferenceOLS一阶差分估计。即对个体固定效应模型中的回归量与被回归量的差分变量构成的模型的参数进行OLS估计。RandomEffects随机效应。当固定效应独立于回归变量时，RandomEffects 使用一个拟差分来有效地估计。 然而，当固定效应和回归变量之间存在依赖关系时，估计量就不再具有一致性。PooledOLS混合最小二乘估计。Pooledols 忽略了固定效应，当效应与回归因子无关时，它是一致的但是低效的。Panelols 比其他的估计方法更通用，可以用于两类效应的模型(例如，实体（即固定）和时间效应)模型中的$\alpha_i$，与回归自变量相关，称之为固定效应模型；与回归自变量不相关，称之为随机效应模型。固定效应将α消掉，随机效应则将其放入误差项，然后探索方差结构。模型规范类似于 statsmodels。 下面例子估计了一个固定效应回归面板上的劳动者的工资，模型中劳动者的工资是经验平方的函数，如果一个男人结婚了虚拟变量设为1，如果这个男人是工会成员同样设定虚拟变量。from linearmodels.panel import PanelOLS from linearmodels.datasets import wage_panel import statsmodels.api as sm data = wage_panel.load() data = data.set_index(['nr','year']) dependent = data.lwage exog = sm.add_constant(data'expersq','married','union') mod = PanelOLS(dependent, exog, entity_effects=True) res = mod.fit(cov_type='unadjusted') res结果： PanelOLS Estimation Summary ================================================================================ Dep. Variable: lwage R-squared: 0.1365 Estimator: PanelOLS R-squared (Between): -0.0674 No. Observations: 4360 R-squared (Within): 0.1365 Date: Wed, Jun 24 2020 R-squared (Overall): 0.0270 Time: 23:03:18 Log-likelihood -1439.0 Cov. Estimator: Unadjusted F-statistic: 200.87 Entities: 545 P-value 0.0000 Avg Obs: 8.0000 Distribution: F(3,3812) Min Obs: 8.0000 Max Obs: 8.0000 F-statistic (robust): 200.87 P-value 0.0000 Time periods: 8 Distribution: F(3,3812) Avg Obs: 545.00 Min Obs: 545.00 Max Obs: 545.00 Parameter Estimates ============================================================================== Parameter Std. Err. T-stat P-value Lower CI Upper CI ------------------------------------------------------------------------------ const 1.3953 0.0123 113.50 0.0000 1.3712 1.4194 expersq 0.0037 0.0002 19.560 0.0000 0.0033 0.0041 married 0.1073 0.0182 5.8992 0.0000 0.0717 0.1430 union 0.0828 0.0198 4.1864 0.0000 0.0440 0.1215 ============================================================================== F-test for Poolability: 9.3360 P-value: 0.0000 Distribution: F(544,3812) Included effects: Entity与statamodel类似，面板模型也可以使用类似 r 的公式来指定。 下面的模型和前面的一模一样。 注意使用特殊变量 EntityEffects 来包含固定效果。mod = PanelOLS.from_formula('lwage ~ 1 + expersq + union + married + EntityEffects',data)数据格式使用面板数据回归模型，需要以下的数据格式二维数据结构：MultiIndex DataFrames，其中外部索引是个体，内部索引是时间，可以使用pandas三维数据结构：3D structures，维度0是变量，维度1是个体索引，维度2是时间索引，还可以使用具有维度(t，n)的2D 数据结构，该数据结构被视为具有维度(1，t，n)的3D 数据结构。 这些3D数据结构可以是 pandas、 NumPy 或 xarray。因为本人习惯用pandas，所以就以pandas为例，其他numpy或者xarry自行去原文档查看。最精确的数据格式是MultiIndex DataFrame，因为只有单个列可以保留面板中的所有类型。下面示例使用培训数据，使用 set_index 命令构造 MultiIndex DataFrame。 个体索引为 fcode，时间索引为year。from linearmodels.datasets import jobtraining data = jobtraining.load() print(data.head()) mi_data = data.set_index(['fcode', 'year']) print(mi_data.head()) employ sales avgsal scrap rework tothrs union grant \ fcode year 410032 1987 100.0 47000000.0 35000.0 NaN NaN 12.0 0 0 1988 131.0 43000000.0 37000.0 NaN NaN 8.0 0 0 1989 123.0 49000000.0 39000.0 NaN NaN 8.0 0 0 410440 1987 12.0 1560000.0 10500.0 NaN NaN 12.0 0 0 1988 13.0 1970000.0 11000.0 NaN NaN 12.0 0 0 d89 d88 ... grant_1 clscrap cgrant clemploy clsales \ fcode year ... 410032 1987 0 0 ... 0 NaN 0 NaN NaN 1988 0 1 ... 0 NaN 0 0.270027 -0.088949 1989 1 0 ... 0 NaN 0 -0.063013 0.130621 410440 1987 0 0 ... 0 NaN 0 NaN NaN 1988 0 1 ... 0 NaN 0 0.080043 0.233347 lavgsal clavgsal cgrant_1 chrsemp clhrsemp fcode year 410032 1987 10.463100 NaN NaN NaN NaN 1988 10.518670 0.055570 0.0 -8.946565 -1.165385 1989 10.571320 0.052644 0.0 0.198597 0.047832 410440 1987 9.259130 NaN NaN NaN NaN 1988 9.305651 0.046520 0.0 0.000000 0.000000 [5 rows x 28 columns]当只引用单个序列时，可以使用 MultiIndex Series 表示形式。from linearmodels import PanelOLS mod = PanelOLS(mi_data.lscrap, mi_data.hrsemp, entity_effects=True) print(mod.fit())案例这些例子涵盖了可用于估计面板模型的模型。 最初的示例都忽略协方差选项，因此使用默认的经典协方差，该协方差适用于 homsedastic （同方差）数据。 备选的协方差选项将在本文档的末尾描述。加载数据这些例子都利用了工资面板数据，参考文献Vella and M. Verbeek (1998), “Whose Wages Do Unions Raise? A Dynamic Model of Unionism and Wage Rate Determination for Young Men,” Journal of Applied Econometrics 13, 163-183.这组数据包括20世纪80年代男性的工资和特征。 个体标识符为nr，标识的时间为year。 这些数据在 Jeffrey Wooldridge 的《计量经济学导论》第14章中得到了广泛的应用。同时，再添加一列year，方便添加时间虚拟变量。from linearmodels.datasets import wage_panel import pandas as pd data = wage_panel.load() year = pd.Categorical(data.year) data = data.set_index(['nr', 'year']) data['year'] = year print(wage_panel.DESCR) print(data.head())nr person identifier year 1980 to 1987 black =1 if black exper labor market experience hisp =1 if Hispanic hours annual hours worked married =1 if married educ years of schooling union =1 if in union lwage log(wage) expersq exper^2 occupation Occupation code基本回归：混合回归仅仅是一个普通的 OLS，方便理解各种面板数据结构。 它作为一个基本模型是很有用的。 在这里，工资的对数被设定为因变量，该模型使用所有的变量和时间虚拟变量。from linearmodels.panel import PooledOLS import statsmodels.api as sm exog_vars = ['black','hisp','exper','expersq','married', 'educ', 'union', 'year'] exog = sm.add_constant(data[exog_vars]) mod = PooledOLS(data.lwage, exog) pooled_res = mod.fit() print(pooled_res)具有不相关效应的参数估计在对面板数据建模时，通常会考虑超出 OLS 有效估计范围的模型。最常见的是误差成分模型，它在标准OLS模型中增加了一个附加项，其中[Math Processing Error(数学处理误差)]影响个体i的所有数值，当[数学处理误差]与[数学处理误差]中的回归因子不相关时，可以使用随机效应模型来有效估计该模型的参数。随机效应RE随机效应模型与混合OLS模型几乎相同，但它考虑到了模型的结构，因此效率更高。随机效应采用了一种quasi-demeaning策略，即减去个体内数值的时间平均数，以考虑到共同的冲击。from linearmodels.panel import RandomEffects mod = RandomEffects(data.lwage, exog) re_res = mod.fit() print(re_res)模型的拟合度相当相似，尽管experience的回归发生了很大变化，其重要性也发生了变化。这在一定程度上是可以解释的，因为加入了年份虚拟变量，会拟合experience的趋势，所以只有横截面的差异才是重要的。随机效应估计器中的准差异取决于一个量，这个量取决于特异性冲击和共同冲击的相对方差。这可以使用variance_decomposition方法来访问。re_res.variance_decompositionEffects 0.106946 Residual 0.123324 Percent due to Effects 0.464438 Name: Variance Decomposition, dtype: float64系数[数学处理误差]决定了demeaning发生的程度。当这个值为1时，RE模型就会还原为pooled模型，因为当effects没有方差时，就会出现这种情况。当面板不平衡时，它会在不同的实体之间变化，但在这个平衡的面板中，所有的值都是一样的。平均数估计BE平均数估计是一种替代性的，通常效率较低的方法，可以用来估计模型参数。它特别简单，因为它首先计算[数学处理误差]和[数学处理误差]的时间平均数，然后用这些平均数进行简单的回归。剔除了年份虚拟变量，因为平均化消除了因年份而产生的差异。 expersq也被剔除，因为它与exper相当共线性。这些结果与前面的模型大致相似。from linearmodels.panel import BetweenOLS exog_vars = ['black','hisp','exper','married', 'educ', 'union'] exog = sm.add_constant(data[exog_vars]) mod = BetweenOLS(data.lwage, exog) be_res = mod.fit() print(be_res)具有相关效应的参数估计FE当效应与回归因子相关时，RE和BE估计量不一致。通常的解决方案是使用PanelOLS中的固定效应。固定效应应用于个体维度时称为实体效应，应用于时间维度时称为时间效应。实体效应通过设置entity_effects=True参数，实体效应被包括在内。这相当于为每个个体加入了虚拟变量。在本面板中，这将增加545个虚拟变量，模型的估计速度将大大降低。PanelOLS实际上并没有使用虚拟变量，而是使用分组demeaning来达到同样的效果。时间不变变量：在使用实体效应时，不能包含时间变量，因为一旦demeaned，这些变量将全部为0。这里也不包括exper，因为一旦纳入实体效应和时间虚拟值，exper是完全共线性的。from linearmodels.panel import PanelOLS exog_vars = ['expersq', 'union', 'married', 'year'] exog = sm.add_constant(data[exog_vars]) mod = PanelOLS(data.lwage, exog, entity_effects=True) fe_res = mod.fit() print(fe_res)时间效应可以使用time_effects=True添加时间效应。这里去掉了时间虚拟变量。请注意，核心系数是相同的。唯一的变化是对集合性的检验统计，因为 "效应 "不同时包括实体和时间，而之前只包括实体。exog_vars = ['expersq','union','married'] exog = sm.add_constant(data[exog_vars]) mod = PanelOLS(data.lwage, exog, entity_effects=True, time_effects=True) fe_te_res = mod.fit() print(fe_te_res)效应VS虚拟变量对于可以一致估计的变量，如通常的大N、固定T面板中的时间效应，没有必要将其作为效应列入。它们可以作为虚拟变量来实施。其他选项当使用PanelOLS和效应时，还有一些额外的选项可以使用。use_lsdv可以用来强制使用虚拟变量来估计模型。这不是必须的，而且在大多数情况下速度会比较慢。auto_df指示PanelOLS确定要进行的自由度调整。与大多数其他估计不同，效果的处理取决于协方差估计器。例如，当使用经典协方差估计器时，效应算作损失的自由度。当使用实体效应和按实体聚类时，则不会。设置auto_df=False将允许实体被计算或不计算，这取决于count_effects的值。一阶差分在可能存在相关性的情况下，一阶差分是固定效应的一种替代方法。在使用一阶差分时，必须排除时间不变的变量。此外，只能包含一个线性时间趋势变量，因为它看起来像一个常数。这个变量将吸收数据中所有的时间趋势，因此对这些变量的解释可能是具有挑战性的。from linearmodels.panel import FirstDifferenceOLS exog_vars = ['exper','expersq', 'union', 'married'] exog = data[exog_vars] mod = FirstDifferenceOLS(data.lwage, exog) fd_res = mod.fit() print(fd_res)比较模型模型结果可以使用compare进行比较。compare接受结果列表、结果字典，其中键被解释为模型名称。from linearmodels.panel import compare print(compare({'BE':be_res,'RE':re_res,'Pooled':pooled_res})) Model Comparison =============================================================== BE RE Pooled --------------------------------------------------------------- Dep. Variable lwage lwage lwage Estimator BetweenOLS RandomEffects PooledOLS No. Observations 545 4360 4360 Cov. Est. Unadjusted Unadjusted Unadjusted R-squared 0.2155 0.1806 0.1893 R-Squared (Within) 0.1141 0.1799 0.1692 R-Squared (Between) 0.2155 0.1853 0.2066 R-Squared (Overall) 0.1686 0.1828 0.1893 F-statistic 24.633 68.409 72.459 P-value (F-stat) 0.0000 0.0000 0.0000 ===================== ============ =============== ============ const 0.2836 0.0234 0.0921 (1.5897) (0.1546) (1.1761) black -0.1414 -0.1394 -0.1392 (-2.8915) (-2.9054) (-5.9049) hisp 0.0100 0.0217 0.0160 (0.2355) (0.5078) (0.7703) exper 0.0278 0.1058 0.0672 (2.4538) (6.8706) (4.9095) married 0.1416 0.0638 0.1083 (3.4346) (3.8035) (6.8997) educ 0.0913 0.0919 0.0913 (8.5159) (8.5744) (17.442) union 0.2587 0.1059 0.1825 (5.6214) (5.9289) (10.635) expersq -0.0047 -0.0024 (-6.8623) (-2.9413) year.1981 0.0404 0.0583 (1.6362) (1.9214) year.1982 0.0309 0.0628 (0.9519) (1.8900) year.1983 0.0202 0.0620 (0.4840) (1.6915) year.1984 0.0430 0.0905 (0.8350) (2.2566) year.1985 0.0577 0.1092 (0.9383) (2.5200) year.1986 0.0918 0.1420 (1.2834) (3.0580) year.1987 0.1348 0.1738 (1.6504) (3.5165) --------------------------------------------------------------- T-stats reported in parentheses协方差选项异方差性稳健共方差Heteroskedasticity Robust CovarianceWhite的稳健协方差可以通过设置cov_type='robust'来使用。这个估计器对某些类型的规范问题增加了一些稳健性，但在使用固定效应（实体效应）时不应使用，因为它不再是稳健的。取而代之的是需要一个聚类协方差。exog_vars = ['black','hisp','exper','expersq','married', 'educ', 'union'] exog = sm.add_constant(data[exog_vars]) mod = PooledOLS(data.lwage, exog) robust = mod.fit(cov_type='robust')按实体聚类通常要聚类的变量是实体或实体和时间。可以使用cov_type='clusterered'和附加关键字参数cluster_entity=True和/或cluster_time=True来实现。clust_entity = mod.fit(cov_type='clustered', cluster_entity=True) clust_entity_time = mod.fit(cov_type='clustered', cluster_entity=True, cluster_time=True)OrderedDict 用于保存比较模型的结果。这允许对模型进行命名，也允许指定模型的顺序。一个标准的 dict 将产生有效的随机顺序。from collections import OrderedDict res = OrderedDict() res['Robust'] = robust res['Entity'] = clust_entity res['Entity-Time'] = clust_entity_time print(compare(res))对实体的聚类降低了整个T统计量。这表明每个实体的残差有重要的相关性。按两者进行聚类也会降低t统计量，这表明数据中存在横截面依赖性。注：按实体聚类解决的是跨时间的相关性，按时间聚类控制的是一个时间段内实体之间的相关性。其他聚类其他聚类可以用输入的聚类直接传递整数组（1列或2列，或1-d数组）。本例按职业进行聚类，由于只有9组，而且通常聚类标准误差的理论要求聚类的数量要大，所以可能不是一个可靠的变量。clust_entity = mod.fit(cov_type='clustered', clusters=data.occupation) print(data.occupation.value_counts()) print(clust_entity)使用公式指定模型所有的模型都可以使用公式来指定，和statsmodels中的公式类似。单变量回归的基础公式语法为y ~ 1 + x用于BetweenOLS、PooledOLS和RandomEffects的公式是完全标准的，与statsmodels相同。FirstDifferenceOLS几乎完全相同，但需要注意的是，该模型不能包含截距。实现效应（实体、时间或其他）的PanelOLS对公式有一个小的扩展，它允许将实体效应或时间效应（或两者）指定为公式的一部分。虽然不能使用公式接口指定其他效应，但在使用公式时，这些效应可以作为一个可选的参数包含在内。加载数据使用公式时，需要一个MultiIndex pandas数据框架，其中的指数是实体时间。这里用statsmodels提供的 "企业投资的决定因素 "中的Grunfeld数据来说明公式的使用。这个数据集包含了企业投资、市场价值和工厂资本存量的数据。from statsmodels.datasets import grunfeld data = grunfeld.load_pandas().data data = data.set_index(['firm','year']) print(data.head())实体效应实体效果使用特殊命令EntityEffects来指定。默认情况下不包含常数，所以如果需要常数，应该在公式中包含1+。当包含效应时，无论是否包含常数，模型和拟合都是相同的。当一个常数被包含在模型中时，并施加额外的约束条件，即效应数为0，这样就可以使用因变量和回归变量的大平均数来确定常数。from linearmodels import PanelOLS mod = PanelOLS.from_formula('invest ~ value + capital + EntityEffects', data=data) print(mod.fit()) mod = PanelOLS.from_formula('invest ~ 1 + value + capital + EntityEffects', data=data) print(mod.fit())实体和时间效应时间效应同样可以使用TimeEffect纳入。在许多模型中，时间效应可以被一致估计，因此它们可以等价地包含在使用分类变量的回归变量集中。mod = PanelOLS.from_formula('invest ~ 1 + value + capital + EntityEffects + TimeEffects', data=data) print(mod.fit())平均数估计其他的面板模型都是直截了当的，为了完整起见，将其包括在内。from linearmodels import BetweenOLS, FirstDifferenceOLS, PooledOLS mod = BetweenOLS.from_formula('invest ~ 1 + value + capital', data=data) print(mod.fit())一阶差分第一个差分模型绝对不能包含一个常数，因为这在差分后是无法确定的。mod = FirstDifferenceOLS.from_formula('invest ~ value + capital', data=data) print(mod.fit())混合回归混合OLS估计是PanelOLS在没有效应时的一种特殊情况。它与statsmodels（或WLS）中的OLS实际上是相同的，但为了完整起见，将其包括在内。mod = PooledOLS.from_formula('invest ~ 1 + value + capital', data=data) print(mod.fit())
前言这其实是我们一次课程作业，以上证50ETF期权为例说明波动率微笑现象。按习惯我先上网搜了一下看有没有前辈写过这样的代码，毕竟重复造轮子不好嘛。没想到真的有，原文链接：https://www.jianshu.com/p/e73f538859df。但是这份代码有个问题，就是需要自己手动搜集数据，而且输出的数据不是标准的DataFrame。趁着做作业的机会，我借鉴并改写了作者的代码，主要实现了以下改进：使用plotly作图，生成可交互式图像。利用tushare自动拉取数据，计算某日所有50ETF期权的隐含波动率，同时可以生成标准的DataFrame数据。借鉴这篇文章的思路，利用看涨看跌期权平价关系计算隐含分红率可生成波动率微笑图和波动率曲面图。详细代码见：https://github.com/caly5144/shu-s-project/tree/master/options理论从期权定价模型本身来说，公式中的波动率指的是未来的波动率数据，这使历史波动率始终存在着较大的缺陷。为了回避这一缺陷，一些学者将目光转向隐含波动率。隐含波动率具体计算思路是：不论用何种理论的理论模型对市场上交易的期权用模型确定其合理的期权价值，一般需要用到6个标准参数：执行价格（K）、到期时间（t）、标的资产价格（S）、利率（rf）、标的资产息票率（q）和波动率。前五个参数可以从市场中和期权合约中获取，只有波动率是未知的参数。因此，将期权市场价格以及除波动率之外的5个参数代入期权定价公式后推导出波动率，由此计算出的波动率称为隐含波动率。因此，隐含波动率是根据期权的市场价格反推出来的波动率，也是市场自身对未来波动率的预期值，是市场参与者通过期权交易的买卖报价对未来波动率达成的共识。上述说明来自蒋祥林《金融衍生工具》数据的话，直接使用tushare即可，这是一个提供各类金融数据的API，我以前也做过介绍，见：https://www.sitstars.com/archives/37/参数tushare函数tushare参数说明执行价格（K）opt_basicexercise_price 到期时间（t）opt_basicdelist_date需要进行处理标的资产价格（S）fund_dailyclose上证50ETF代码为510050.SH利率（rf）shibor1w其实用7日银行间国债回购利率更好，但是tushare不提供，只能用一周shibor利率代替期权最新价格（c）opt_dailysettle 分红率用看涨看跌期权平价关系求出：$$ c\left(t\right)+Ke^{-r\left(T-t\right)}=p\left(t\right)+e^{-q\left(T-t\right)}S\left(t\right) $$将上述参数代入Black-Scholes定价公式，即可求出相应隐含波动率$$ d_1=\frac{\ln{\left(S\left(t\right)/K\right)}+\left(r-q+\sigma^2/2\right)\left(T-t\right)}{\sigma\sqrt{T-t}} $$$$ d_2=\frac{\ln{\left(S\left(t\right)/K\right)}+\left(r-q-\sigma^2/2\right)\left(T-t\right)}{\sigma\sqrt{T-t}}=d_1-\sigma\sqrt{T-t} $$看涨期权：$$ c\left(t\right)=S\left(t\right)e^{-q\left(T-t\right)}N\left(d_1\right)-Ke^{-r\left(T-t\right)}N\left(d_2\right) $$看跌期权：$$ p=Ke^{-r\left(T-t\right)}N\left(-d_2\right)-Se^{-q\left(T-t\right)}N\left(-d_1\right) $$代码实现导入库import pandas as pd import datetime import tushare as ts import numpy as np from math import log,sqrt,exp from scipy import stats import plotly.graph_objects as go import plotly import plotly.express as px pro = ts.pro_api() plotly.offline.init_notebook_mode(connected=True)提取数据并清洗def extra_data(date): # 提取数据 # 提取50ETF合约基础信息 df_basic = pro.opt_basic(exchange='SSE', fields='ts_code,name,call_put,exercise_price,list_date,delist_date') df_basic = df_basic.loc[df_basic['name'].str.contains('50ETF')] df_basic = df_basic[(df_basic.list_date<=date)&(df_basic.delist_date>date)] # 提取当天市场上交易的期权合约 df_basic = df_basic.drop(['name','list_date'],axis=1) df_basic['date'] = date # 提取日线行情数据 df_cal = pro.trade_cal(exchange='SSE', cal_date=date, fields = 'cal_date,is_open,pretrade_date') if df_cal.iloc[0, 1] == 0: date = df_cal.iloc[0, 2] # 判断当天是否为交易日，若否则选择前一个交易日 opt_list = df_basic['ts_code'].tolist() # 获取50ETF期权合约列表 df_daily = pro.opt_daily(trade_date=date,exchange = 'SSE',fields='ts_code,trade_date,settle') df_daily = df_daily[df_daily['ts_code'].isin(opt_list)] # 提取50etf指数数据 df_50etf = pro.fund_daily(ts_code='510050.SH', trade_date = date,fields = 'close') s = df_50etf.iloc[0, 0] # 提取无风险利率数据（用一周shibor利率表示） df_shibor = pro.shibor(date = date,fields = '1w') rf = df_shibor.iloc[0,0]/100 # 数据合并 df = pd.merge(df_basic,df_daily,how='left',on=['ts_code']) df['s'] = s df['r'] = rf df = df.rename(columns={'exercise_price':'k', 'settle':'c'}) #print(df) return df def data_clear(df): # 数据清洗 def days(df): # 计算期权到期时间 start_date = datetime.datetime.strptime(df.date,"%Y%m%d") end_date = datetime.datetime.strptime(df.delist_date,"%Y%m%d") delta = end_date - start_date return int(delta.days)/365 def iq(df): # 计算隐含分红率 #q = -log((df.settle+df.exercise_price*exp(-df.interest*df.delta)-df.settle_p)/(df.s0))/df.delta q = -log((df.c+df.k*exp(-df.r*df.t)-df.c_p)/(df.s))/df.t return q df['t'] = df.apply(days,axis = 1) df = df.drop(['delist_date','date'],axis = 1) # 计算隐含分红率 df_c = df[df['call_put']=='C'] df_p = df[df['call_put']=='P'] df_p = df_p.rename(columns={'c':'c_p','ts_code':'ts_code_p', 'call_put':'call_put_p'}) df = pd.merge(df_c,df_p,how='left',on=['trade_date','k','t','r','s']) df['q'] = df.apply(iq,axis = 1) c_list = [x for x in range(8)]+[11] df_c = df.iloc[:,c_list] df_p = df'ts_code_p','trade_date','c_p','call_put_p','k','t','r','s','q' df_p = df_p.rename(columns={'c_p':'c','ts_code_p':'ts_code', 'call_put_p':'call_put'}) df_c = df_c.append(df_p) #print(df_c) return df_c 计算隐含波动率因为公式中有正态分布存在，所以我们没办法直接求出解，只能用二分法求出数值解。#根据BS公式计算期权价值 def bsm_value(s,k,t,r,sigma,q,option_type): d1 = ( log( s/k ) + ( r -q + 0.5*sigma**2 )*t )/( sigma*sqrt(t) ) d2 = ( log( s/k ) + ( r -q - 0.5*sigma**2 )*t )/( sigma*sqrt(t) ) if option_type.lower() == 'c': value = (s*exp(-q*t)*stats.norm.cdf( d1) - k*exp( -r*t )*stats.norm.cdf( d2)) else: value = k * exp(-r * t) * stats.norm.cdf(-d2) - s*exp(-q*t) * stats.norm.cdf(-d1) return value #二分法求隐含波动率 def bsm_imp_vol(s,k,t,r,c,q,option_type): c_est = 0 # 期权价格估计值 top = 1 #波动率上限 floor = 0 #波动率下限 sigma = ( floor + top )/2 #波动率初始值 count = 0 # 计数器 while abs( c - c_est ) > 0.000001: c_est = bsm_value(s,k,t,r,sigma,q,option_type) #根据价格判断波动率是被低估还是高估，并对波动率做修正 count += 1 if count > 100: # 时间价值为0的期权是算不出隐含波动率的，因此迭代到一定次数就不再迭代了 sigma = 0 break if c - c_est > 0: #f(x)>0 floor = sigma sigma = ( sigma + top )/2 else: top = sigma sigma = ( sigma + floor )/2 return sigma def cal_iv(df): # 计算主程序 option_list = df.ts_code.tolist() df = df.set_index('ts_code') alist = [] for option in option_list: s = df.loc[option,'s'] k = df.loc[option,'k'] t = df.loc[option,'t'] r = df.loc[option,'r'] c = df.loc[option,'c'] q = df.loc[option,'q'] option_type = df.loc[option,'call_put'] sigma = bsm_imp_vol(s,k,t,r,c,q,option_type) alist.append(sigma) df['iv'] = alist return df 多项式拟合填充缺失值由于部分期限的50etf期权品种较少，所以使用多项式拟合的方法补全这部分数据，由波动率微笑可知，这是一个二次函数，因此令degree = 2此外，用plotly做3D曲面图时，需要数据格式如下： 因此，用pandas的pivot_table函数进行整理。def data_pivot(df): # 数据透视 df = df.reset_index() option_type = 'C' # 具有相同执行价格、相同剩余到期时间的看涨看跌期权隐含波动率相等，因此算一个就够了 df = df[df['call_put']==option_type] df = df.drop(['ts_code','trade_date','c','s','r','call_put','q'],axis=1) df['t'] = df['t']*365 df['t'] = df['t'].astype(int) df = df.pivot_table(index=["k"],columns=["t"],values=["iv"]) df.columns = df.columns.droplevel(0) df.index.name = None df = df.reset_index() df = df.rename(columns={'index':'k'}) return df def fitting(df): # 多项式拟合 col_list = df.columns for i in range(df.shape[1]-1): x_col = col_list[0] y_col = col_list[i+1] df1 = df.dropna(subset=[y_col]) x = df1.iloc[:,0] y = df1.iloc[:,i+1] degree = 2 weights = np.polyfit(x, y, degree) model = np.poly1d(weights) predict = np.poly1d(model) x_given_list = df[pd.isnull(df[y_col]) == True][x_col].tolist() # 所有空值对应的k组成列表 for x_given in x_given_list: y_predict = predict(x_given) df.loc[df[x_col]==x_given, y_col] = y_predict return df作图def im_surface(): # 波动率曲面作图 df = plot_df() df = fitting(df) #df.to_excel('iv_fitting.xlsx') df = df.set_index('k') y = np.array(df.index) x = np.array(df.columns) fig = go.Figure(data=[go.Surface(z=df.values, x=x, y=y)]) fig.update_layout(scene = dict( xaxis_title='剩余期限', yaxis_title='执行价格', zaxis_title='隐含波动率'), width=1400, margin=dict(r=20, b=10, l=10, t=10)) #fig.write_image("fig1.jpg") plotly.offline.plot(fig) def smile_plot(): # 波动率微笑作图 df = plot_df() df = df.set_index('k') df = df.stack().reset_index() df.columns = ['k', 'days', 'iv'] fig = px.line(df, x="k", y="iv", color="days",line_shape="spline") plotly.offline.plot(fig)综合def main(): date = '20200515' # 可任意更改日期 df = extra_data(date) # 提取数据 df = data_clear(df) # 数据清洗 df = cal_iv(df) # 计算隐含波动率 data_pivot(df) # 数据透视表 smile_plot() # 波动率微笑 im_surface() # 波动率曲面 if __name__ == '__main__': main()效果
这是用python做计量经济学实证分析系列第二篇。原文地址：https://python.quantecon.org/mle.html[scode type="green"]本系列文章对原文进行了一定的修改，主要是作图工具从matplotlib变为plotly.express，因为后者更为简洁明了，如果你不习惯plotly，可以看原文的代码[/scode][scode type="yellow"]本系列文章需要你对计量经济学有一定的了解，并且熟悉python、pandas、numpy[/scode]引入在上一节课中，我们使用线性回归来估计因变量和解释变量之间的关系。但如果我们的模型并不是线性关系呢？一种广泛使用的替代方法是最大似然估计，它包括指定一类由未知参数索引的分布，然后使用数据来确定这些参数值。相对于线性回归的好处是，它使变量之间的概率关系具有更大的灵活性。这里我们通过复制丹尼尔·特雷斯曼（Daniel Treisman，2016）的论文《俄罗斯的亿万富翁》（Russia's Billionals）来说明最大似然估计，该论文将一个国家的亿万富翁数量与其经济特征联系起来。本文的结论是，俄罗斯的亿万富翁数量高于市场规模和税率等经济因素的预测。文献地址：http://pubs.aeaweb.org/doi/pdfplus/10.1257/aer.p20161068导入库：import numpy as np from numpy import exp from scipy.special import factorial import pandas as pd from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.api import Poisson from scipy import stats from scipy.stats import norm from statsmodels.iolib.summary2 import summary_col import plotly.express as px import plotly plotly.offline.init_notebook_mode(connected=True)开始简介最大似然估计的第一步是选择数据的概率分布。更准确地说，我们需要假设参数服从什么分布类。如正态分布或Gamma分布。每一类都是由有限个参数索引的分布族。如正态分布是一类均值为$\mu \in(-\infty, \infty)$，标准差$\sigma \in(0, \infty)$的正态分布。我们将通过固定参数的分布让数据挑选出特定元素。由此产生的参数估计称为最大似然估计。亿万富翁Treisman对估算不同国家的亿万富翁数量很感兴趣。亿万富翁的数量是整数值。因此我们考虑只取非负整数的分布。（这就是为什么最小二乘回归不是解决当前问题的最佳工具，因为线性回归中的因变量不限于整数值）一个整数分布是Poisson分布，其概率质量函数（probability mass function，pmf）为：$${f(y)}=\frac{\mu^{y}}{y !} e^{-\mu}, \quad y=0,1,2, \ldots, \infty$$对于u的不同值，我们可以绘制y上的Poisson分布，如下所示poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ) df_fin = pd.DataFrame() for μ in [1, 5, 10]: df = pd.DataFrame() distribution = [poisson_pmf(y_i, μ) for y_i in range(0, 25)] df['y'] = list(range(0, 25)) df['f'] = distribution df['μ'] = μ df_fin = df_fin.append(df) fig = px.line(df_fin, x="y", y="f",color="μ",labels = {'f':'f(y/μ)'}) fig.update_traces(mode='markers+lines') plotly.offline.plot(fig) 注意，随着y的平均值的增加，Poisson分布开始类似于正态分布。让我们来看看我们在这堂课中要处理的数据的分布情况。Treisman的主要数据来源是《福布斯》对亿万富翁及其估计净值的排名。原数据集在这：https://www.aeaweb.org/articles?id=10.1257/aer.p20161068df = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/mle/fp.dta') print(df.head()) 使用柱状图，我们可以看到2008年每个国家亿万富翁人数的分布，(numbil0)（出于作图目的舍弃了美国的数据）numbil0_2008 = df[(df['year'] == 2008) & (df['country'] != 'United States')] fig = px.histogram(numbil0_2008, x="numbil0",labels={'numbil0':'Number of billionaires in 2008'}) plotly.offline.plot(fig) 从直方图来看，泊松假设并非不合理（尽管u值很低，且存在一些异常值）条件分布在Treisman的论文中，因变量——i国亿万富翁的数量$y_i$——它被模拟为人均国内生产总值、人口规模以及加入关贸总协定和世贸组织的年数的函数。因此，yi的分布需要以解释变量$x_i$的向量为条件。标准的poisson回归模型如下：$$\begin{array}{c}f\left(y_{i} | \mathbf{x}_{i}\right)=\frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}} ; \quad y_{i}=0,1,2, \ldots, \infty \qquad(1)\\ \text { where } \mu_{i}=\exp \left(\mathbf{x}_{i}^{\prime} \beta\right)=\exp \left(\beta_{0}+\beta_{1} x_{i 1}+\ldots+\beta_{k} x_{i k}\right)\end{array}$$为了说明yi的分布依赖于xi，让我们做一个简单的拟合。使用上面的poisson_pmf函数，并代入$\beta$和$x_i$的任意值。poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ) df_fin = pd.DataFrame() # Define a parameter vector with estimates β = np.array([0.26, 0.18, 0.25, -0.1, -0.22]) # Create some observations X datasets = [np.array([0, 1, 1, 1, 2]), np.array([2, 3, 2, 4, 0]), np.array([3, 4, 5, 3, 2]), np.array([6, 5, 4, 4, 7])] for X in datasets: μ = round(exp(X @ β),3) distribution = [poisson_pmf(y_i, μ) for y_i in range(0, 25)] df = pd.DataFrame() df['y'] = list(range(0, 25)) df['f'] = distribution df['μ'] = μ df_fin = df_fin.append(df) fig = px.line(df_fin, x="y", y="f",color="μ", labels = {'f':r'$f(y \mid x_i; \beta)$',,'y':'$y \mid x_i$'}) fig.update_traces(mode='markers+lines') plotly.offline.plot(fig,include_mathjax='cdn') 我们可以看到yi的分布是以xi为条件的；（ui不再是常数）最大似然估计在我们的亿万富翁模型中，条件分布包含4个我们需要估计的参数（k=4）。我们将整个参数向量标记为$\beta$，其中$$\beta=\left[\begin{array}{l}\beta_{0} \\ \beta_{1} \\ \beta_{2} \\ \beta_{3}\end{array}\right]$$为了用极大似然估计法估计模型，我们希望使我们的估计$\hat\beta$是真参数$\beta$的可能性最大。首先，我们需要构造似然函数（$\mathcal{L}(\beta)$），它类似于联合概率密度函数(probability density function，pdf)。假设我们有一些数据$y_{i}=\left\{y_{1}, y_{2}\right\} \text { and } y_{i} \sim f\left(y_{i}\right)$。如果y1和y2是独立的，则这些数据的联合pmf是$f\left(y_{1}, y_{2}\right)=f\left(y_{1}\right) \cdot f\left(y_{2}\right)$。如果$y_i$服从Poisson分布，$\lambda=7$，我们可以可视化联合pmf如下：import matplotlib.pyplot as plt def plot_joint_poisson(μ=7, y_n=20): poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ) yi_values = np.arange(0, y_n, 1) # Create coordinate points of X and Y X, Y = np.meshgrid(yi_values, yi_values) # Multiply distributions together Z = poisson_pmf(X, μ) * poisson_pmf(Y, μ) fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') ax.plot_surface(X, Y, Z.T, cmap='terrain', alpha=0.6) ax.scatter(X, Y, Z.T, color='black', alpha=0.5, linewidths=1) ax.set(xlabel='$y_1$', ylabel='$y_2$') ax.set_zlabel('$f(y_1, y_2)$', labelpad=10) plt.show() plot_joint_poisson(μ=7, y_n=20)不会用plotly做类似的图，只能照抄源码用matplotlib.pyplot了。 类似地，我们数据的联合pmf（作为条件Poisson分布）可以写成$$f\left(y_{1}, y_{2}, \ldots, y_{n} | \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n} ; \beta\right)=\prod_{i=1}^{n} \frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}}$$yi取决于$x_i$和参数$\beta$的值。似然函数与联合PMF相同，但将参数$\beta$作为随机变量对待，并给出观测值（$y_i$，$x_i$）。$$\begin{aligned} \mathcal{L}\left(\beta | y_{1}, y_{2}, \ldots, y_{n} ; \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n}\right) &=\prod_{i=1}^{n} \frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}} \\ &=f\left(y_{1}, y_{2}, \ldots, y_{n} | \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n} ; \beta\right) \end{aligned}$$既然我们有了似然函数，我们想找到产生最大似然值的$\beta$$$\max _{\beta} \mathcal{L}(\beta)$$log-likelihood通常更容易最大化（考虑区分${f(x)}=x \exp (x) \text { vs. } {f(x)}=\log (x)+x$）假设取对数是单调递增变换，则似然函数的极大值也是对数似然函数的极大值。在我们的例子中，log-likelihood是：$$\begin{aligned} \log \mathcal{L}(\beta)=& \log \left(f\left(y_{1} ; \beta\right) \cdot f\left(y_{2} ; \beta\right) \cdot \ldots \cdot f\left(y_{n} ; \beta\right)\right) \\ &=\sum_{i=1}^{n} \log f\left(y_{i} ; \beta\right) \\ &=\sum_{i=1}^{n} \log \left(\frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}}\right) \\ &=\sum_{i=1}^{n} y_{i} \log \mu_{i}-\sum_{i=1}^{n} \mu_{i}-\sum_{i=1}^{n} \log y ! \end{aligned}$$然而，上述问题不存在解析解-要找到最大似然估计，我们需要使用数值方法。数值方法许多分布没有很好的解析解，因此需要数值方法来求解参数估计。其中一种数值方法是牛顿-拉斐逊算法。我们的目标是找到最大似然估计$\hat\beta$。在$\hat\beta$处，log-likelihoo函数的一阶导数等于0。我们假设：$$\log \mathcal{L}(\beta)=-(\beta-10)^{2}-10$$β = np.linspace(1, 20) logL = -(β - 10) ** 2 - 10 dlogL = -2 * β + 20 df = pd.DataFrame({'logL':logL,'dlogL':dlogL,'β':β}) df_long=pd.melt(df, id_vars=['β'], value_vars=['logL', 'dlogL']) fig = px.line(df_long, x="β", y="value", facet_row="variable",line_shape="spline") fig['layout']['yaxis']['title']['text']=r'$\frac{dlog \mathcal{L(\beta)}}{d \beta}$' fig['layout']['yaxis2']['title']['text']=r'$log \mathcal{L(\beta)}$' fig['layout']['annotations'][0]['text'] = '' fig['layout']['annotations'][1]['text'] = '' fig.update_yaxes(matches=None) plotly.offline.plot(fig,include_mathjax='cdn') 上图显示了最大似然值（顶部图表）出现在$\frac{d \log \mathcal{L}(\beta)}{d \beta}=0$时。因此，当$\beta=10$时，似然值最大。我们还可以通过检查二阶导数（底部曲线图的斜率）是否为负值来确保该值是最大值（而不是最小值）。Newton-Raphson算法找到一个一阶导数为0的点。为了使用该算法，我们在最大值处进行初始猜测（OLS参数估计可能是合理的猜测），然后1.使用更新规则迭代算法$$\beta_{(k+1)}=\beta_{(k)}-H^{-1}\left(\beta_{(k)}\right) G\left(\beta_{(k)}\right)$$其中：$$\begin{aligned} G\left(\beta_{(k)}\right) &=\frac{d \log \mathcal{L}\left(\beta_{(k)}\right)}{d \beta_{(k)}} \\ H\left(\beta_{(k)}\right) &=\frac{d^{2} \log \mathcal{L}\left(\beta_{(k)}\right)}{d \beta_{(k)} d \beta_{(k)}^{\prime}} \end{aligned}$$2.检查是否$\beta_{(k+1)}-\beta_{(k)}<t o l$如果为true，则停止迭代并设置$\hat{\beta}=\beta_{(k+1)}$如果为false，则更新$\beta_{(k+1)}$从更新方程可以看出，只有当$G\left(\beta_{(k)}\right)=0$即一阶导数等于0时，$\beta_{(k+1)}=\beta_{(k)}$，（实际上，当差异低于一个可容忍的阈值时，我们就停止迭代）让我们开始实现Newton-Raphson算法。首先，我们将创建一个名为poissonrege的类，这样我们就可以很容易地重新计算每次迭代的对数似然、梯度和Hessian值class PoissonRegression: def __init__(self, y, X, β): self.X = X self.n, self.k = X.shape # Reshape y as a n_by_1 column vector self.y = y.reshape(self.n,1) # Reshape β as a k_by_1 column vector self.β = β.reshape(self.k,1) def μ(self): return np.exp(self.X @ self.β) def logL(self): y = self.y μ = self.μ() return np.sum(y * np.log(μ) - μ - np.log(factorial(y))) def G(self): X = self.X y = self.y μ = self.μ() return X.T @ (y - μ) def H(self): X = self.X μ = self.μ() return -(X.T @ (μ * X))我们的函数newton_raphson将取一个PoissonRegression 对象，该对象具有对参数向量$\beta_0$的初始猜测。该算法根据更新规则更新参数向量，并在新的参数估计下重新计算梯度矩阵和Hessian矩阵。迭代将在下列任一情况下结束：参数和更新参数之间的差异低于可容忍水平。已达到最大迭代次数（意味着未达到收敛）因此，我们可以了解算法运行时的情况，在每次迭代时添加一个选项display=True来打印值。def newton_raphson(model, tol=1e-3, max_iter=1000, display=True): i = 0 error = 100 # Initial error value # Print header of output if display: header = f'{"Iteration_k":<13}{"Log-likelihood":<16}{"θ":<60}' print(header) print("-" * len(header)) # While loop runs while any value in error is greater # than the tolerance until max iterations are reached while np.any(error > tol) and i < max_iter: H, G = model.H(), model.G() β_new = model.β - (np.linalg.inv(H) @ G) error = β_new - model.β model.β = β_new # Print iterations if display: β_list = [f'{t:.3}' for t in list(model.β.flatten())] update = f'{i:<13}{model.logL():<16.8}{β_list}' print(update) i += 1 print(f'Number of iterations: {i}') print(f'β_hat = {model.β.flatten()}') # Return a flat array for β (instead of a k_by_1 column vector) return model.β.flatten()让我们用一个X包含5个观测值和3个变量的小数据集来验证我们的算法。X = np.array([[1, 2, 5], [1, 1, 3], [1, 4, 2], [1, 5, 2], [1, 3, 1]]) y = np.array([1, 0, 1, 1, 0]) # Take a guess at initial βs init_β = np.array([0.1, 0.1, 0.1]) # Create an object with Poisson model values poi = PoissonRegression(y, X, β=init_β) # Use newton_raphson to find the MLE β_hat = newton_raphson(poi, display=True) 由于这是一个观测较少的简单模型，6次迭代后就收敛了。可以看到，随着每次迭代，log似然值都会增加。记住，我们的目标是最大化对数似然函数，这是算法一直致力于实现的。另外，注意随着迭代次数的增加，$\log \mathcal{L}\left(\beta_{(k)}\right)$的增加值变小了。这是因为当我们达到最大值时，梯度接近0，因此我们的更新方程中的分子变得更小。梯度向量在$\hat\beta$处应该接近于0poi.G() >>> [[-3.952e-07] [-1.001e-06] [-7.731e-07]]迭代过程可以在下图中可视化，其中最大值出现在$\beta=10$处logL = lambda x: -(x - 10) ** 2 - 10 def find_tangent(β, a=0.01): y1 = logL(β) y2 = logL(β+a) x = np.array([[β, 1], [β+a, 1]]) m, c = np.linalg.lstsq(x, np.array([y1, y2]), rcond=None)[0] return m, c β = np.linspace(2, 18) fig, ax = plt.subplots(figsize=(12, 8)) ax.plot(β, logL(β), lw=2, c='black') for β in [7, 8.5, 9.5, 10]: β_line = np.linspace(β-2, β+2) m, c = find_tangent(β) y = m * β_line + c ax.plot(β_line, y, '-', c='purple', alpha=0.8) ax.text(β+2.05, y[-1], f'$G({β}) = {abs(m):.0f}$', fontsize=12) ax.vlines(β, -24, logL(β), linestyles='--', alpha=0.5) ax.hlines(logL(β), 6, β, linestyles='--', alpha=0.5) ax.set(ylim=(-24, -4), xlim=(6, 13)) ax.set_xlabel(r'$\beta$', fontsize=15) ax.set_ylabel(r'$log \mathcal{L(\beta)}$',rotation=0,labelpad=25,fontsize=15) ax.grid(alpha=0.3) plt.show() 请注意，我们对Newton-Raphson算法的实现是相当基本的—有关更健壮的实现，请参见scipy.optimize。使用statsmodels进行最大似然估计知道了背后的原理，我们可以将MLE应用到一个有趣的应用程序中。我们将在statsmodels中使用Poisson回归模型来获得更丰富的输出，包括标准误差、测试值等等。statsmodels使用与上面相同的算法来寻找最大似然估计。在我们开始之前，让我们用stats模型重新估计我们的简单模型，以确认我们得到了相同的系数和对数似然值。X = np.array([[1, 2, 5], [1, 1, 3], [1, 4, 2], [1, 5, 2], [1, 3, 1]]) y = np.array([1, 0, 1, 1, 0]) stats_poisson = Poisson(y, X).fit() print(stats_poisson.summary()) 现在，让我们复制丹尼尔·特雷斯曼（Daniel Treisman）在报告中提到的俄罗斯亿万富翁论文的结果。Treisman从估计等式（1）开始，式中：yi为亿万富翁数$x_{i1}$是人均GDP的对数$x_{i2}$是人口对数$x_{i3}$是加入GATT和世贸组织时长（代表进入国际市场）本文只考虑2008年的估计。我们将像这样设置我们的估计变量df = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/mle/fp.dta') df = df[df['year'] == 2008] df['const'] = 1 # Variable sets reg1 = ['const', 'lngdppc', 'lnpop', 'gattwto08'] reg2 = ['const', 'lngdppc', 'lnpop','gattwto08', 'lnmcap08', 'rintr', 'topint08'] reg3 = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08','rintr', 'topint08', 'nrrents', 'roflaw']然后利用statsmodels中的Poisson函数对模型进行拟合。在作者的论文中，我们将使用健壮的标准错误# Specify model poisson_reg = sm.Poisson(df'numbil0', df[reg1],missing='drop').fit(cov_type='HC0') print(poisson_reg.summary()) 成功！算法在9次迭代中收敛。我们的结果表明，人均国内生产总值（GDP）和加入关税及贸易总协定（GATT）的年限与一个国家的亿万富翁数量正相关，正如预期的那样。让我们来估计一下作者的更全功能的模型，并将它们显示在一个表中regs = [reg1, reg2, reg3] reg_names = ['Model 1', 'Model 2', 'Model 3'] info_dict = {'Pseudo R-squared': lambda x: f"{x.prsquared:.2f}", 'No. observations': lambda x: f"{int(x.nobs):d}"} regressor_order = ['const','lngdppc','lnpop','gattwto08','lnmcap08','rintr', 'topint08','nrrents','roflaw'] results = [] for reg in regs: result = sm.Poisson(df'numbil0', df[reg], missing='drop').fit(cov_type='HC0',maxiter=100, disp=0) results.append(result) results_table = summary_col(results=results,float_format='%0.3f',stars=True, model_names=reg_names,info_dict=info_dict, regressor_order=regressor_order) results_table.add_title('Table 1 - Explaining the Number of Billionaires in 2008') print(results_table) 研究结果表明，亿万富翁的数量与人均GDP、人口规模、股市市值呈正相关，与最高边际所得税率呈负相关。为了按国家分析我们的结果，我们可以绘制出预测值和实际值之间的差异，然后从最高值到最低值排序，并绘制出前15个国家data = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08','rintr', 'topint08', 'nrrents', 'roflaw','numbil0','country'] results_df = df[data].dropna() # Use last model (model 3) results_df['prediction'] = results[-1].predict() # Calculate difference results_df['difference'] = results_df['numbil0'] - results_df['prediction'] # Sort in descending order results_df = results_df.sort_values('difference', ascending=False)[:15] # Plot the first 15 data points fig = px.bar(results_df, x="country", y="difference") plotly.offline.plot(fig,include_mathjax='cdn') 正如我们所看到的，俄罗斯的亿万富翁人数远远超过了该模型的预测（比预期多出约50人）Treisman利用这一实证结果讨论了俄罗斯亿万富翁过剩的可能原因，包括俄罗斯财富的来源、政治气候以及苏联解体后的私有化历史。总结在本课程中，我们使用最大似然估计来估计泊松模型的参数。statsmodels包含其他内置的似然模型，如Probit和Logit。可以在这里找到示例：https://www.statsmodels.org/dev/examples/notebooks/generated/generic_mle.html
这是用python做计量经济学实证分析系列第一篇。原文地址：https://python.quantecon.org/ols.html[scode type="green"]本系列文章对原文进行了一定的修改，主要是作图工具从matplotlib变为plotly.express，因为后者更为简洁明了，如果你不习惯plotly，可以看原文的代码[/scode][scode type="yellow"]本系列文章需要你对计量经济学有一定的了解，并且熟悉python、pandas、numpy[/scode]引入文献：Acemoglu D, Johnson S, Robinson J A. The colonial origins of comparative development: An empirical investigation[J]. American economic review, 2001, 91(5): 1369-1401.地址：https://economics.mit.edu/files/4123本文强调了制度在经济发展中的重要性。其主要贡献是利用殖民者死亡率作为制度差异外生变化的来源。这样的变化是必要的，以确定制度是否导致了更快的经济增长，而不是相反。安装相应库：pip install linearmodels导入相应库import numpy as np import pandas as pd import plotly.express as px import plotly import statsmodels.api as sm from statsmodels.iolib.summary2 import summary_col from linearmodels.iv import IV2SLS plotly.offline.init_notebook_mode(connected=True)简单的线性回归Acemoglu研究的是制度差异是否导致了经济产出的差异。但是，我们如何衡量制度和产出呢？在论文中，经济产出以经汇率调整后的1995年人均GDP的对数代表。制度差异以Political Risk Services Group构建的1985-95年期间的平均征收保护指数（protection against expropriation on average）代表。数据可以在Daron Acemoglu的个人网站上下载。https://economics.mit.edu/faculty/acemoglu/data/ajr2001我们使用read_stata() 读取.dta格式文件。df1 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable1.dta') print(df1.head()) 对logGDP和平均征收保护指数（avexpr）做散点图fig = px.scatter(df1,x='avexpr',y='logpgp95') plotly.offline.plot(fig) 上图显示了征收保护指数与人均GDP的对数之间相当强的正相关关系。具体来说，如果更高的征收保护是衡量制度质量的一个指标，那么更好的制度似乎与更好的经济成果（人均国内生产总值）呈正相关。选择一个线性模型来描述这种关系似乎是合理的假设。我们可以把模型写成$$\operatorname{logpgp} 95_{i}=\beta_{0}+\beta_{1} \text { avexpr }_{i}+u_{i}$$其中，$\beta_{0}$是趋势线在y轴上的截距$\beta_{1}$是线性趋势线的斜率，表示征收保护指数对log人均GDP的边际影响u是一个随机误差项（由于模型中未包含的因素，观测值与线性趋势有偏差）绘制趋势线# 删除空值 df1_subset = df1.dropna(subset=['logpgp95', 'avexpr']) # Use only 'base sample' for plotting purposes df1_subset = df1_subset[df1_subset['baseco'] == 1] fig = px.scatter(df1_subset, x="avexpr", y="logpgp95", trendline="ols", title={'text': "Figure 2: OLS relationship between expropriation risk and income",'y':0.9,'x':0.5}, text="shortnam") fig.update_traces(textposition='top center') fig.update_xaxes(title_text='Average Expropriation Risk 1985-95') fig.update_yaxes(title_text='Log GDP per capita, PPP, 1995') plotly.offline.plot(fig) 最常用的估计$\beta$的方法是最小二乘法（OLS），顾名思义，OLS模型是通过寻找使残差平方和最小化的参数来求解的。即：$$\min _{\hat{\beta}} \sum_{i=1}^{N} \widehat{u}_{i}^{2}$$为了求解常数项$\beta_0$，我们需要在数据集中添加一列1（即用$\beta_0x_i,x_i=1$来替换$\beta_0$）df1['const'] = 1现在我们可以使用statsmodels中的OLS函数构建我们的模型。reg1 = sm.OLS(endog=df1['logpgp95'], exog=df1'const', 'avexpr', \ missing='drop') results = reg1.fit() # 获取参数的估计值 print(results.summary()) # 打印回归结果 从上图中，我们可以看到截距$\hat \beta_0=4.63$斜率$\hat\beta_1 = 0.53$$\hat\beta_1$为正，意味着制度质量对经济成果有正向作用。$\hat\beta_1$的 p-value为0，意味着上述结论在统计学上是显著的（使用p<0.05为拒绝条件）R平方值为0.611表明，人均国内生产总值的对数61%左右的变化可以由征收保护指数解释。使用我们的参数估计，我们现在可以将我们的估计关系写为$$\hat{logpgp} 95_{i}=4.63+0.53 \text { avexpr }_{i}$$我们可以用这个方程来预测征收保护指数指数对的人均GDP水平的影响。例如，对于一个指数值为7.07（数据集的平均值）的国家，我们发现其1995年人均对数国内生产总值的预测水平为8.38。mean_expr = np.mean(df1_subset['avexpr']) predicted_logpdp95 = 4.63 + 0.53 * 7.07获得这个结果的一个更简单（更准确）的方法是使用predict()方法，设置常数项constant = 1 and $avexpr_i= mean\_expr$results.predict(exog=[1, mean_expr]) # array([8.09156367])同样地，我们可以使用predict()预测所有值作图并与实际值比较。# Drop missing observations from whole sample df1_plot = df1.dropna(subset=['logpgp95', 'avexpr']) df1_plot['predict'] = results.predict() df_long=pd.melt(df1_plot, id_vars=['avexpr'], value_vars=['logpgp95', 'predict']) fig = px.scatter(df_long, x='avexpr', y='value', color='variable', title={'text': "OLS predict",'y':0.9,'x':0.5}, labels={'logpgp95':'observed'},symbol="variable") fig.update_yaxes(title_text='logpgp95') plotly.offline.plot(fig) 线性回归的扩展-多元回归到目前为止，我们只考虑了影响经济表现的制度因素——几乎可以肯定的是，还有许多其他影响GDP的因素没有包括在我们的模型中。剔除影响logpgp95的其他变量将导致遗漏变量偏差（omitted variable bias），产生有偏差和不一致的参数估计。我们可以通过加入其他可能影响logpgp95i的因素，将我们的二元回归模型扩展到多元回归模型。考虑其他因素例如：气候对经济表现的影响，用纬度来代表气候经济表现和制度（如文化、历史等）的差异；用大陆虚拟变量（dummies）来控制让我们使用maketable2.dta中的数据估计本文中考虑的一些扩展模型df2 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable2.dta') df2['const'] = 1 # 创建参数列表 X1 = ['const', 'avexpr'] X2 = ['const', 'avexpr', 'lat_abst'] X3 = ['const', 'avexpr', 'lat_abst', 'asia', 'africa', 'other'] # OLS回归 reg1 = sm.OLS(df2['logpgp95'], df2[X1], missing='drop').fit() reg2 = sm.OLS(df2['logpgp95'], df2[X2], missing='drop').fit() reg3 = sm.OLS(df2['logpgp95'], df2[X3], missing='drop').fit()现在我们已经拟合了模型，我们将使用summary_col在一个表中显示结果（模型编号与本文中的编号相对应）info_dict={'R-squared': lambda x:f"{x.rsquared:.2f}", 'No.observations': lambda x: f"{int(x.nobs):d}"} # f是去掉转义字符，nobs就是总观测数量 results_table = summary_col(results=[reg1,reg2,reg3],float_format='%0.2f', stars = True,model_names=['Model 1','Model 3','Model 4'], info_dict=info_dict, regressor_order=['const','avexpr','lat_abst','asia','africa']) results_table.add_title('Table 2 - OLS Regressions') print(results_table) 内生性-两阶段最小二乘OLS模型可能存在内生性问题，从而导致模型估计的偏差和不一致。也就是说，制度与经济成果之间可能存在双向关系：较富裕的国家可能负担得起或更喜欢更好的制度影响收入的变量也可能与制度差异相关指数的构建可能有失偏颇：研究者可能倾向于研究收入较高、拥有更好制度的国家为了处理内生性，我们可以使用两阶段最小二乘（2SLS）回归，这是OLS回归的扩展。此方法需要用以下变量替换内生变量avexpr：与avexpr相关与误差项无关（即不应直接影响因变量，否则由于omitted variable bias，他会与$u_i$相关；）新的回归变量被称为工具变量，其目的是消除我们的制度差异代理中的内生性。原文献的主要贡献是利用殖民者的死亡率来衡量制度差异。他们假设，殖民者的高死亡率导致建立了更具extractive in nature的制度（对征用的保护更少，这些制度至今仍然存在）。使用散点图，我们可以看到征收保护与殖民者死亡率呈负相关，符合作者的假设并满足有效工具变量的第一个条件。df1_subset2 = df1.dropna(subset=['logpgp95', 'avexpr']) # Use only 'base sample' for plotting purposes df1_subset2 = df1.dropna(subset=['logem4', 'avexpr']) fig = px.scatter(df1_subset2, x="logem4", y="avexpr", trendline="ols", title={'text': "Figure 3: First-stage relationship between settler mortality and expropriation risk",'y':0.9,'x':0.5}, text="shortnam") fig.update_traces(textposition='top center') fig.update_xaxes(title_text='Log of Settler Mortality') fig.update_yaxes(title_text='Average Expropriation Risk 1985-95') plotly.offline.plot(fig) 如果17至19世纪的殖民者死亡率对当前的国内生产总值有直接影响（除了通过制度产生的间接影响外），第二个条件可能无法满足例如，定居者死亡率可能与一个国家当前的疾病环境有关，这可能影响当前的经济表现。原文献认为这不太可能，因为：大多数定居者死亡是由于疟疾和黄热病，对当地居民的影响有限。例如，非洲或印度当地人民的疾病负担似乎没有高于平均水平，这是殖民前这些地区相对较高的人口密度造成的。这似乎是一个有效的工具，我们可以使用2SLS回归获得一致和无偏的参数估计。第一阶段第一阶段的因变量是内生变量（$avexpr_i$）。工具是我们模型中所有外生变量的集合（而不仅仅是我们替换的变量）。以模型为例，我们的工具只是一个常数和殖民者死亡率对数$\log e m 4_{i}$，因此，我们将第一阶段回归估计为$$\text { avexpr }_{i}=\delta_{0}+\delta_{1} \log e m 4_{i}+v_{i}$$我们需要估计这个方程的数据位于maketable4.dta中（仅用baseco = 1的数据进行估算）df4 = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/ols/maketable4.dta') df4 = df4[df4['baseco'] == 1] # Add a constant variable df4['const'] = 1 # Fit the first stage regression and print summary results_fs = sm.OLS(df4['avexpr'],df4'const', 'logem4',missing='drop').fit() print(results_fs.summary()) 二阶段我们需要使用predict()检索$\text {avexpr}_{i}$的预测值。然后，我们将原始线性模型中的内生变量$\text {avexpr}_{i}$替换为预测值$\widehat {avexpr}_{i}$。因此，我们的第二阶段回归为：$$\operatorname{logpgp} 95_{i}=\beta_{0}+\beta_{1} a \widehat{v e x p} r_{i}+u_{i}$$df4['predicted_avexpr'] = results_fs.predict() results_ss = sm.OLS(df4['logpgp95'],df4'const', 'predicted_avexpr').fit() print(results_ss.summary()) 第二阶段的回归结果为我们提供了一个制度对经济结果的影响的无偏且一致的估计。结果表明，与OLS结果相比，两者之间存在更强的正相关关系。请注意，虽然我们的参数估计是正确的，但标准误不是，因此，不建议“手动”（使用OLS分阶段）计算2SLS。我们可以使用linearmodels包一步正确地进行2SLS回归。注意，当使用IV2SLS时，外生变量和工具变量在函数参数中是分开的iv = IV2SLS(dependent=df4['logpgp95'],exog=df4['const'],endog=df4['avexpr'], instruments=df4['logem4']).fit(cov_type='unadjusted') print(iv.summary) 鉴于我们现在有一致和无偏见的估计，我们可以从我们估计的模型中推断，制度差异（源自殖民时期建立的制度）有助于解释当今各国收入水平的差异。原文献使用0.94的边际效应来计算智利和尼日利亚之间的指数差异（即制度质量）意味着收入差异高达7倍，强调了制度在经济发展中的重要性。