Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈点点 (中国人民大学)
邮箱:chendiandian@ruc.edu.cn
编者按:本文主要摘译自下文,特此致谢!
Source:Jinjing L I, Zyphur M, Sugihara G. EDM: Stata module to implement empirical dynamic modeling[J]. 2019. -Link-
目录
如何研究不满足经典回归假设的复杂动态系统?在这些系统中,准自然实验或基于模型的回归往往无法适用,因果识别更是一项困难的任务。处理现实世界的复杂性需要能够描述和测试非线性动态系统因果关系的工具。
为了解决现有方法的局限性,研究人员开发了一种新方法,即经验动态建模 (empirical dynamic modeling, EDM)。这种来自自然科学的新方法允许:
Stata 中的 edm
命令可以实现这一方法。该命令为时间序列和面板数据提供了三种经验动态建模 (EDM) 方法:
接下来,我们将使用模拟数据来说明这些方法。
EDM 的逻辑基于这样一个事实,即产生数据的动态系统可以通过重构其基础系统来建模 (Takens 1981)。建议观看以下视频以直观了解 EDM 的思想:
考虑一个具有 D 个时变变量的系统,例如一个沿着 GDP、失业和失业率变化的国民经济。随着时间的推移,D 个变量绘制了系统状态变化的轨迹图。下图是 D=3 情况的示意图 (呈现出 “蝴蝶” 吸引子)。随着国民经济发展,D 个变量的轨迹将在 D 维状态空间中形成 D 维流形
如果吸引子流形
使用滞后来重构流形是一种状态空间重构的 “延迟嵌入” 或 “滞后坐标嵌入” 方法,其中
为了将
另一半数据形成 “预测” 或测试 / 验证集,其中包含落在
每个邻居 i 的权重由目标到该邻居的欧几里德距离和距离衰减参数
其中,
预测质量通过预测集的目标的未来实现值与投影单形的加权平均值的相关性
为了推断维度为 D 的基础系统,在不同的 E 值下评估
理想情况下,一个系统可以用少于 10 个因素 (即,少于 10 个维度) 来描述,这样当 E<10 时预测最准确。在这种情况下,系统可被视为低维和确定性系统,其预测精度很高 (例如,
与单纯形投影一样,默认情况下,S - 映射使用 50/50 的数据分割为 E 长度嵌入向量的训练集和预测集。该训练集表示一个重构流形
其中,系数向量
当权重函数中的
再次,使用
与单纯形投影一样,S 映射也可以以多变量方式进行。在多变量情况下,当
收敛交叉映射基于这样一个事实,即如果 X 是 Y 的确定性驱动因素 (X→Y),则 Y 的状态必须包含有助于恢复或 “交叉映射” X 状态的信息。该方法是单纯形投影的一个扩展,即使用一个变量重构吸引子流形 M,并使用它来预测不同的变量。如果变量共享同一个吸引子流形,则可以使用重构的流形进行预测。
具体来说,如果变量 X 和 Y 参与同一个动态系统流形 M,则重构的流形
这是违反直觉的,因为在典型的时间序列或面板数据方法中,原因是用来解释或预测结果,而不是相反。然而,在 CCM 中,结果变量 Y 交叉映射或 “xmaps” 原因变量 X,通过阴影流形
CCM 中的术语 “收敛” 描述了评估因果关系的标准。该术语反映了这样一个事实,即如果存在 X→Y 因果关系,则预测精度 (例如,X 和
下面我们通过一个常见的动态系统来说明 EDM 命令的使用。两个变量 x 和 y 都依赖于各自的过去值,并通过方程的最后一项相互耦合。本例中,x 设置为仅由其自身的过去值确定,而 y 由 y 和 x 的过去值确定。
下图绘制了
EDM 分析的第一步是确定系统的维度。维度可以理解为重建基础流形 explore
子命令的单纯形投影完成的,使用 e()
选项可以指定维数范围。
使用 explore
子命令,将数据随机分成两半,其中一半用作训练数据集以构造阴影流形
在本例中,我们使用单纯形投影探索 2 到 10 之间的所有维度,以确定最佳 E。对于大型数据集,可以选择 E=20 或更高。由于训练数据集是随机选择的,我们将该方法复制了 50 次,通过选项 rep(50)
来估计 50 次运行的平均性能。命令和输出如下所示。ci()
选项可用于报告从多次复制中得出的 edm
命令也可与 jackknife pre fix
结合使用,以实现附加控制。通常选择具有最高
. edm explore y, e(2/10) rep(50)
Replication progress (50 in total)
.................................................. 50
Empirical Dynamic Modelling
Univariate mapping with y and its lag values
----------------------------------------------------------------------
--------- rho --------- --------- MAE ---------
Actual E theta Mean Std. Dev. Mean Std. Dev.
----------------------------------------------------------------------
2 1 .97999 .0059908 .031708 .004511
3 1 .97083 .0078036 .039437 .0044651
4 1 .95058 .013928 .049484 .0060995
5 1 .92494 .023275 .059684 .0078145
6 1 .89361 .032127 .069651 .0094557
7 1 .8451 .046727 .082331 .010843
8 1 .79302 .052859 .096854 .010539
9 1 .75721 .058412 .10613 .010663
10 1 .72319 .063342 .11468 .010283
----------------------------------------------------------------------
Note: Results from 50 runs
Note: Number of neighbours (k) is set to E+1
Note: Random 50/50 split for training and validation data
结果表明,随着 E 的增加,e(explore_result)
矩阵的内容绘制结果,如下图所示。
如果流形 M 上的轨迹相对于系统当前状态不变,则存在线性;而如果系统演化依赖于状态,则存在非线性。通过从单纯形投影中选择最佳 E 后,接下来我们可以估计一个回归模型,该模型使用距离衰减参数
当选项 theta(0)
中
explore
和 xmap
子命令都支持 S 映射,而不是局部预测的单纯形投影。下面给出了一个例子来分析前一个例子中变量 y 的非线性。在这种情况下,我们探索 0 和 5 之间可能的 k()
选项中指定一个负数来包含本地预测的所有观测值。这允许在低 E 或
. edm explore y, e(2) algorithm(smap) theta(0(0.01)5) k(-1)
Empirical Dynamic Modelling
Univariate mapping with y and its lag values
--------------------------------------------------------------------
Actual E theta rho MAE
--------------------------------------------------------------------
2 0 .67288 .13505
2 .01 .67831 .13412
2 .02 .68365 .1332
2 .03 .6889 .13228
(output omitted)
2 4.94 .99378 .019535
2 4.95 .99379 .019527
2 4.96 .99379 .01952
2 4.97 .9938 .019513
2 4.98 .99381 .019505
2 4.99 .99382 .019498
2 5 .99382 .01949
--------------------------------------------------------------------
Note: Number of neighbours (k) is adjusted to 74
Note: Random 50/50 split for training and validation data
上面命令的结果也可以使用返回矩阵 e(explore_result)
以图形方式绘制,如下图所示。
在前面的步骤中,分析表明本例的最佳 E 为 2,与数据生成过程一致。现在,我们使用 xmap
子命令,使用之前观察到的 e=2 值,推导变量 x 和 y 之间的双变量 (整体) 因果效应。虽然在大多数情况下,我们建议将变量标准化,使其处于相同的尺度,但在这个假设的示例中,标准化对结果影响不大。
CCM 使用一个变量的重构流形 (潜在结果) 预测另一个变量 (潜在原因) 的值来评估变量之间的偶然联系。本质上,这是一种双变量单纯形投影,使用一个变量形成训练集,同时对另一个变量进行预测。使用预测值与观测值,计算
值得注意的是,单纯形投影结果可能表明不同的变量应使用不同的 E 值。如果是这种情况,选择用于重建流形的 E 可与 direction(oneway)
选项一起使用,以对每个变量运行单独的 CCM 程序。在这种情况下,应该采用原因变量的 E,因为原因变量是由结果预测的。例如,如果对 y 而言 E=2,对 x 而言 E=4,则可以运行以下 CCM:edm xmap x y,e (2) direction (oneway)
,用于 edm xmap y x,e (4) direction (oneway)
用于
. edm xmap x y, e(2)
Empirical Dynamic Modelling
Convergent Cross-mapping result for variables x and y
---------------------------------------------------------------
Mapping Library size rho MAE
---------------------------------------------------------------
y ~ y|M(x) 150 .23019 .19673
x ~ x|M(y) 150 .69682 .13714
---------------------------------------------------------------
Note: The embedding dimension E is 2
默认情况下,在不为 library()
指定值的情况下,整个数据集将用作训练库。换句话说,来自一个变量的所有可能信息都用于预测另一个变量。这通常是 CCM 的第一步,因为事先不容易知道最优 L 是什么。此外,如果最大 L 的 CCM 没有显示出有意义的可预测性,那么收敛的问题就没有意义了。然后,可以在评估整个结果时手动设定最大 L。具体而言,如果存在因果关系 library()
选项指定了库大小的范围,replicate()
选项允许重复估计。当指定的 L 小于最大可用库长时,将为重建流形选择一个随机子样本,并且可以选择 replicate()
选项。
下面的示例估计了库长为 5 到 150 时的 xmap
过程,然后在较小的库长下使用复制选项来评估不确定性。
. edm xmap x y, e(2) rep(10) library(5/150)
Replication progress (20 in total)
....................
Empirical Dynamic Modelling
Convergent Cross-mapping result for variables x and y
---------------------------------------------------------------
Mapping Lib size Mean rho Std. Dev.
---------------------------------------------------------------
y ~ y|M(x) 5 -.0090982 .12325
y ~ y|M(x) 6 -.004514 .11706
y ~ y|M(x) 7 .018604 .10854
y ~ y|M(x) 8 .01794 .10815
y ~ y|M(x) 9 .033454 .095494
y ~ y|M(x) 10 .037815 .097865
y ~ y|M(x) 11 .030339 .10197
y ~ y|M(x) 12 .040244 .10324
y ~ y|M(x) 13 .049327 .098762
(output omitted)
x ~ x|M(y) 142 .68351 .010373
x ~ x|M(y) 143 .68683 .0093171
x ~ x|M(y) 144 .68719 .0092867
x ~ x|M(y) 145 .6898 .011021
x ~ x|M(y) 146 .68971 .0098815
x ~ x|M(y) 147 .6911 .0073416
x ~ x|M(y) 148 .69069 .0061557
x ~ x|M(y) 149 .69297 .0054071
x ~ x|M(y) 150 .69682 0
---------------------------------------------------------------
Note: Results from 10 replications
Note: The embedding dimension E is 2
详细的结果存储在两个返回矩阵 e(xmap_1)
和 e(xmap_2)
中,可以使用 Stata 的 svmat
命令将其转换为变量,以便于操作。e(xmap_1)
将变量 x 视为变量 y 的结果 (即 e(xmap_2)
将变量 y 视为变量 x (即
还应注意,当两种预测相似且处于预测性能都较好时,结果不排除双向因果关系的可能性。在这种情况下,预测性能的相对差异可用于确定更主要的因果方向,而不排除双向因果关系的推断。当从两个方向的预测性能都较差时,则变量可能是因果独立的。当两个变量的可预测性都很高,但没有观察到收敛,那么有可能存在一个外生的第三个变量可在驱动这两个变量。
/* edm paper replication do file
Jinjing Li, University of Canberra
Michael Zyphur, University of Melbourne
*/
clear all
set more off
set scheme sj
/* Example with a Synthetic Dataset */
/* Create a dynamic system */
set obs 500
gen t = _n
tsset t
gen x = 0.2 if _n==1
gen y = 0.3 if _n==1
gen z = 0.1 if _n==1
local r_x 3.79
local r_y 3.79
local r_z 3.77
local beta_xy = 0.0
local beta_yx=0.2
local beta_zy = 0.5
local tau = 1
drawnorm u1 u2
qui {
forvalues i=2/`=_N' {
replace x=l.x *(`r_x' *(1-l.x)-`beta_xy'*l.y) in `i'
replace y=l.y *(`r_y' *(1-l.y)-`beta_yx'*l`tau'.x) in `i'
replace z=l.z *(`r_z' *(1-l.z)-`beta_zy'*l`tau'.y) in `i'
}
}
keep in 300/450
set seed 12345678
/* Typical Output */
edm explore x
edm xmap x y
/* Plot the system */
twoway (connected x t)(connected y t) in 1/40, ytitle(value) xtitle(t) ///
legend(position(7) ring(0))
graph export plot.pdf, replace
/* rho-E plot */
edm explore y, e(2/10) rep(50)
mat r= e(explore_result)
svmat r, names(col)
twoway (scatter c3 c1)(lpoly c3 c1),xtitle("E") ytitle("{it:{&rho}}") ///
legend(order(1 "{it:{&rho}}" 2 "local polynomial smoothing") ///
col(1) position(8) ring(0))
graph export rho-e.pdf, replace
drop c*
/* rho-theta plot */
edm explore y, e(2) algorithm(smap) theta(0(0.01)5) k(-1)
mat r= e(explore_result)
svmat r, names(col)
twoway (line c3 c2) , legend(order(1 "{it:{&rho}}") position(5) ring(0)) ///
xtitle("{it:{&theta}}") ytitle("{it:{&rho}}") ///
title("{it:{&rho}}-{it:{&theta}} of variable y")
graph export rho-theta.pdf, replace
drop c*
/* ccm */
edm xmap x y, e(2)
/* ccm-convergent plot */
edm xmap x y, e(2) rep(10) library(5/150)
mat c1= e(xmap_1)
mat c2= e(xmap_2)
svmat c1, names(xy)
svmat c2, names(yx)
label variable xy3 "y|M(x)"
label variable yx3 "x|M(y)"
twoway (scatter xy3 xy2, mfcolor(%30) mlcolor(%30)) ///
(scatter yx3 yx2, mfcolor(%30) mlcolor(%30)) ///
(lpoly xy3 xy2)(lpoly yx3 yx2), xtitle(L) ytitle("{it:{&rho}}") legend(col(2))
graph export rho-L.pdf, replace
Note:产生如下推文列表的 Stata 命令为:
lianxh 动态 系统 edm, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh