Shu's Garden - 计量学习
type: channel
link: https://www.sitstars.com/category/%E8%AE%A1%E9%87%8F%E5%AD%A6%E4%B9%A0/
items: 【译】用python做计量之面板数据模型
items: 使用python计算期权隐含波动率并作图
items: 【译】使用python做计量之最大似然估计
items: 【译】使用python做计量之线性回归

【译】用python做计量之面板数据模型
type: item
date: 2020-06-25T01:27:00+08:00
description: 原文：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 Parame…

使用python计算期权隐含波动率并作图
type: item
date: 2020-05-21T23:58:00+08:00
description: 前言这其实是我们一次课程作业，以上证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合约基础…

【译】使用python做计量之最大似然估计
type: item
date: 2020-05-12T11:46:48+08:00
description: 这是用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.ae…

【译】使用python做计量之线性回归
type: item
date: 2020-05-12T11:38:42+08:00
description: 这是用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"…