R语言金融状态空间模型.pdf

1、4.状 态 空 间 模 型 /state space model 李 智 Q群 ： 261821669 July 15, 20160.1 状 状 状 态 态 态 空 空 空 间 间 间 模 模 模 型 型 型 这 里 是 Tsay book第 11章 的 内 容 。 让 y 1 代 表 黄 金 的 市 场 报 价 ， 让 x 1 代 表 黄 金 的 真 实 价 值 ， 我 们 可 以 有 下 面 的 线 性 关 系 ： y 1 =a 1 x 1 +v 1 : (1) 首 先 a 1 代 表 斜 率 。 其 次 这 个 线 性 关 系 中 省 略 了 截 距 ， 因 为 是 否 放 入 截 距

2、是 要 具 体 问 题 具 体 分 析 的 ， 这 里 先 省 略 它 为 的 是 大 家 看 得 清 楚 。 再 次 v 1 是 一 个 normal error。 我 们 只 能 观 测 到 黄 金 的 市 场 报 价 ， 黄 金 的 真 实 价 值 是 观 测 不 到 的 ， 所 以 黄 金 的 真 实 价 值 是 所 谓 的 隐 藏 状 态 the hidden state。 考 虑 多 元 模 型 ， 让 y 1 代 表 黄 金 的 市 场 报 价 ， y 2 代 表 白 银 的 市 场 报 价 ， y 3 代 表 铂 金 的 市 场 报 价 。 让 x 1 代 表 黄 金 的 真

3、实 价 值 ， x 2 代 表 白 银 的 真 实 价 值 ， x 3 代 表 铂 金 的 真 实 价 值 。 我 们 可 以 建 立 两 种 相 似但 不 同 的 线 性 关 系 ： y 1 =a 1 x 1 +v 1 (2) y 2 =a 2 x 2 +v 2 y 3 =a 3 x 3 +v 3 或 ， y 1 =a 1;1 x 1 +a 1;2 x 2 +a 1;3 x 3 +v 1 (3) y 2 =a 2;1 x 2 +a 2;2 x 2 +a 2;3 x 3 +v 2 y 3 =a 3;1 x 3 +a 3;2 x 2 +a 3;3 x 3 +v 3 : 我 们 倾 向 于 使 用

4、 Equation(3)作 为 模 型 ， 因 为 它 能 更 充 分 地 利 用 信 息 和 更 精 确 地 模 型 这 个 线 性 系 统 。 所 以 我 们 就 需 要 把 Equation(3)用 向 量 和 矩 阵 的 方 式 重 新 书 写出 来 ， 成 为 y =Ax +v; y是 向 量 fy 1 ;y 2 ;y 3 g， x是 向 量 fx 1 ;x 2 ;x 3 g， A是 系 数 矩 阵 fa ij g， v是 向 量 fv 1 ;v 2 ;v 3 g normal errors。 黄 金 ， 白 银 ， 铂 金 是 一个 portfolio。 加 上 时 间 t， y

5、t 向 量 是 这 个 投 资 组 合 在 某 一 时 间 点 的 市 场 报 价 。 x t 向 量 是 这 个 投 资 组 合 在 这 个 时 间 点 上 的 真 实 价 值 。 我 们 假 设 ： 第 一 ， 市 场 报 价 是 由 真 实 价 值 决 定 的 ， 市 场 报 价 和 真 实 价 值 之 间 有 线 性 关 系 。 第 二 ， 现 在 的 真 实 价 值 是 并 且 只 是 由 前 一个 时 间 节 点 的 真 实 价 值 决 定 的 ， 所 以 现 在 的 真 实 价 值 和 前 一个 时 间 节 点 的 真 实 价 值 也 存 在 着 线 性 关 系 。 基 于 这

6、两个 假 设 ， 我 们 就 有 了一个 最 简 单 的 状 态 空 间 模 型 ： y t =A t x t +v t measurement/observation equation (4) x t = x t 1 +w t transition/state equation (5) 1其 中 ， 是 现 在 的 真 实 价 值 和 前 一 时 间 点 真 实 价 值 线 性 关 系 的 系 数 矩 阵 ， 它 没 有 t作 下 标 ， 说 明 假 设 它 是 不 随 时 间 变化 的 。 还 有 w t 是 向 量 fw 1 ;w 2 ;w 3 g另 外 的 一些 normal erro

7、rs。 黄 金 ， 白 银 ， 铂 金 的 真 实 价 值 是 不 能 被 观 测 到 的 ， 它 们 被 称 为不 可 见 的 状 态 the hidden states。 在 上 面 的 两 条 基 本 假 设 里 ， 我 们 说 现 在 的 状 态 x t 只 是 由 前 一个 状 态 x t 1 所 决 定 ， 这 完 全 符 合 AR(1)时 间 序 列 模 型 的 定 义 ， 而 且 这 同 样 也 符 合 马 科 夫 过 程 (Markov Process)的 定 义 。 在 状 态 空 间 模 型 里 (state space models)， 能 够 作 为 数 据 观 测

8、到 的 只 是 y t ， 例 如 前 面 提 到 的 市 场 报 价 ， y t 通 常 被 称 作 measurement。 所 以 ， 状 态 空 间 模 型 的目的 就 是 通 过 观 测 measurement， 建 立 模 型 ， 再 估 计 模 型 里 的 系 数 矩 阵 ， 不 论 有 多 少 个 ， 或 者 多复 杂 ， 从 而 估 计 这 个 模 型 里 每 个 时 间 上 的 状 态 向 量 x t 。 值 得 注 意 的 是 ， 现 在 估 计 的 状 态 向 量 x t 还 要 输 入 到 下一个 计 算 中 ， 计 算 下一个 状 态 向 量 x t+1 ， 这 被

9、 称 作 反 馈 (feedback)。 一 般 来 讲 ， 那 些 类 似 于 A t , 的 系 数 矩 阵 是 通 过 最 大 似 值 估 计 (MLE)计 算 的 ， 状 态 向 量 x t 的 估 计 是 通 过 Kalman lter 计 算 的 。 计 算 过 去 已 经 发 生 过 的 状 态 叫 做 smoothing， 计 算 现 在 的 状 态 叫 做 ltering ， 计 算 未来 的 状 态 叫 做 forecasting。 当 然 ， Kalman lter 就 是 计 算 ltering 。 我 真 的 不 愿意 进 入 Kalman lter 的 细 节 ，

10、因 为 符 合 我 们 前 面 那 两个 基 本 假 设 的 状 态 空 间 模 型 有 很 多 也 复 杂 的 多 ， 它 们 的 Kalman lter 都 必 须 具 体 问 题 具 体 分 析 。 0.2 dlm 在 R语 言 里 ， 空 间 状 态 模 型 主 要 是 用 dlm包 做 的 (dynamic linear models)。 可 以 做 state space ARMAX， state space ARIMA， state space Stochastic Volatility， 这 些 都 在 书中 有 详 细 记 载 ， 有 需 要 的 可 以 参 考 。 Line

11、ar regres- sion model也 可 以 写 成 空 间 状 态 的 形式 ， 而 且 估 计 出 来 的 截 距 和 斜 率 是 随 着 时 间 的 变化 而 变化 的 time varying parameters(tvp)。 首 先 ， state space linear regression的 写 法 如 下 ： y t = t + t x t +v t ; v t N(0; 2 t ) (6) t = t 1 +w ;t w ;t N(0; 2 ) (7) t = t 1 +w ;t w ;t N(0; 2 ) (8) 上 面 的 t ; t 是 the hidden

12、states， 可 以 把 它 们 合 并 成 一个 状 态 向 量 t = f t ; t g 0 ， 我 们 的 向 量 都 是 竖 着 的 ， 所 以 要 加 个 撇 ， 竖 过 来 。 v t 是 属 于 y t 的 normal error， w ;t 是 属 于 t 的 normal error， w ;t 是 属 于 t 的 normal error。 下 面 是 state space linear regression的 矩 阵 运 算 形式 ： 2y t = (1 x t ) t +v t (9) t t = 1 0 0 1 t 1 t 1 + w ;t w ;t (10)

13、 我 们 看 到 了 前 面 提 到 过 的 A t 矩 阵 (1 x t )， 和 矩 阵 1 0 0 1 。 下一 步 ， 我 们 把 Equation(9,10)输 入 到 dlm环 境 里 。 然 后 ， 调 取 两个 股 票 的 数 据 序 列 XOM和 JPM， 生 成 一个 空 间 状 态 模 型 ， 再 做 ltering ， 再 做 smoothing， 再 做 forecasting。 XOM是 y t ， JMP是 x t 。 图 Figure (1)输 出 的 是 随 时 间 变化 的 斜 率 和 截 距 ， 可 以 看 到 是 非 常 有 规 律 的 。 图 Figu

14、re (2)输 出 的 是 预 测 的 未来 的 斜 率 和 截 距 。 require(quantmod) require(dlm) require(PerformanceAnalytics) symbols - c(“XOM“,“JPM“) # portfolio symbolData - new.env() getSymbols(symbols, env = symbolData, src = “yahoo“, from=“2012-01-01“, to=“2012-12-31“) ClosePrices - do.call(cbind,eapply(symbolData, Cl) xo

15、m - ClosePrices,2 jpm - ClosePrices,1 buildTVP - function(parm, x.mat) parm - exp(parm) return( dlmModReg(X=x.mat, dV=parm1, dW=c(parm2, parm3) ) start.vals = c(0,0,0) names(start.vals) = c(“lns2v“, “lns2a“, “lns2b“) TVP.mle = dlmMLE(y=xom, parm=start.vals, x.mat=jpm, build=buildTVP, hessian=T) se2

16、- sqrt(exp(TVP.mle$par) names(se2) = c(“sv“, “sa“, “sb“) # fitted ss model 3Figure 1:空 间 状 态 Linear Regression随 时 间 变化 的 斜 率 和 截 距 4Figure 2:空 间 状 态 Linear Regression斜 率 和 截 距 的 预 测 5TVP.dlm - buildTVP(TVP.mle$par, jpm) # filtering TVP.f - dlmFilter(xom, TVP.dlm) # smoothing TVP.s - dlmSmooth(TVP.f)

17、 # extract smoothed states - intercept and slope coefs alpha.s = xts(TVP.s$s-1,1,drop=FALSE, as.Date(rownames(TVP.s$s-1,) beta.s = xts(TVP.s$s-1,2,drop=FALSE, as.Date(rownames(TVP.s$s-1,) colnames(alpha.s) = “alpha“ colnames(beta.s) = “beta“ # extract std errors - dlmSvd2var gives list of MSE matric

18、es mse.list = dlmSvd2var(TVP.s$U.S, TVP.s$D.S) se.mat = t(sapply(mse.list, FUN=function(x) sqrt(diag(x) se.xts = xts(se.mat-1, , index(beta.s) colnames(se.xts) = c(“alpha“, “beta“) a.u = alpha.s + 1.96*se.xts,“alpha“ a.l = alpha.s - 1.96*se.xts, “alpha“ b.u = beta.s + 1.96*se.xts,“beta“ b.l = beta.s

19、 - 1.96*se.xts, “beta“ layout(1:2) chart.TimeSeries(cbind(alpha.s, a.l, a.u), main=“Smoothed estimates of alpha“, colorset=c(1,2,2), lty=c(1,2,2),ylab=expression(alpha),xlab=“) chart.TimeSeries(cbind(beta.s, b.l, b.u), main=“Smoothed estimates of beta“, colorset=c(1,2,2), lty=c(1,2,2),ylab=expressio

20、n(beta),xlab=“) # forecasting using dlmFilter # add 10 missing values to end of sample new.xts = xts(rep(NA, 10), seq.Date(from=end(xom), by=“months“, length.out=11)-1) xom.ext = merge(xom, new.xts),1 TVP.ext.f = dlmFilter(xom.ext, TVP.dlm) # extract h-step ahead forecasts of state vector TVP.ext.f$

21、mas.character(index(new.xts), 60.3 Exponential Smoothing 状 态 空 间 模 型 给 我 们 的 启 示 是 ， 前 面 的 线 性 系 统 Equation(4, 5)是 没 有 解 的 ， MLE只 是 一个 求 解 的 方 法 ， 当 然 还 有 其 他 的 解 法 。 我 们 必 须 要 给 出 一个 空 间 状 态 模 型 交 易 的 例 子 ， 用 一些 R代 码 。 这 里 要 做 的 是 空 间 状 态 Exponential Smoothing,一个 很 好 的 预 测 模 型 ， 可 以 用 来 帮 助 我 们交 易

22、。 要 调 用 forecast包 ， 使 用 其 中 的 ets()函 数 (error trend seasons)。 这 个 函 数 非 常 好 用 ， 所 有 的 参 数 都 是 自 动 估 计 的 ， 而 且 使 用 模 型 的 类 型 也 是 自 动 根 据 AIC BIC法 则 选 取 的 。 如 果 历史 数 据 里 有 trend， 他 就 会 自 动 在 模 型 里 加 入 trend组 件 ， 如 果有 season他 就 会 自 动加 入 season组 件 。 我 们 通 常 称 error trend seasons是 micro structures。 我 们 使

23、 用 的 具 体 策 略 是 ， 使 用 30天 的 历史 数 据 days=30， 预 测 未来 一 天 的 回 报 率 。 如 果 今 天 得 到 的 回 报 率 小 于 预 测 的 明 天 的 回 报 率 ， 我 们 就 在 明 天 做 长 ， 反 而 言 之 ， 在 明 天 做 短 。 这 样 就 需 要 用 一个 函 数 来 决 定 是 否 在 明 天 做 长 还 是 短 。 这 个 函 数 叫 做 regime， 只 有 3行 命 令 。 函 数 的 输 入 叫 做 rolling， 其 实 只 是 30天 的 历史 回 报 率 。 第 一 行 ， 用 ets()函 数 建 模 。

24、 第 二 行 ， 把 模 型 输 入 到 forecast()函 数 里 做 向 前 一 步 的 预 测 。 第 三 行 ， 如 果 今 天 的 回 报 率 小 于 预 测 的 回 报 率 ， 输 出 1， 反 之 输 出 -1。 我 们 需 要 每 天 做 一 个 预 测 ， 然 后 一 天 一 天 的 往 前 推 ， 则 在 R语 言 里 叫 做 rollapply， 也 就 是 使 用 zoo包 ， 把 regime()函 数 植 入 rollapply()函 数 就 可 以了 。 数 列 r存 储 的 是 每 一 天 的 回 报 率 。 把 数 列 r和 regime()函 数 植 入

25、 rollapply， R语 言 就 会 自 动 在 数 列 r上 执 行 regime()函 数 。 如 果 我 们 预 测 明 天 的 回 报 率 会 高 于今 天 的 回 报 率 ， 决 定 做 长 ， 就 先 做 了个 标 记 1。 然 后 在 第 二 天 的 回 报 率 是 .3， 说 明 我 们 当 天 投 资 的 方 向 是 对 的 ， 当 天 投 资 的 回 报 率 是 1 :3 =:3。 再 比 如 ， 如 果 预 测 明 天 回 报 率 会低 于今 天 的 回 报 率 ， 决 定 做 短 ， 标 记 -1。 然 后 在 第 二 天 回 报 率 是 -.3， 我 们 当 天

26、的 投 资 方 向 又 是 正 确 的 ， 当 天 的 投 资 回 报 率 是 1 ( :3) =:3。 函 数 regime()输 出 的 是 投 资 方 向 ， rollapply把 这 些 投 资 方 向 存 入 side数 列 。 我 们 把所 以 相 应 的 第 二 天 的 SP500回 报 率 存 入 next.r数 列 。 然 后 把 它 们 乘 起 来 ， 就 是 每 一 天 的 投 资 回 报 ： trades =side next:r (11) 回 报 率 累 积 的 结 果 是 收 益 prot ， 回 报 率 的 累 积 是 一个 简 单 乘 法 关 系 ， 就 是 银

27、 行 里 计 算 利 息 ： profit = n Y i=1 (trades i + 1) (12) 最 后 我 们 还 需 要 做 两 张 图 Figure: (3)， 上 面 的 是 SP500的 价 格 线 ， 下 面 的 是 使 用 这 个 简 单 策 略 在 这 段 时 间 里 的 收 益 线 。 我 们 使 用 的 数 据 是 2007-01- 01到 2009-01-01 SP500价 格 。 收 益 最 终 达 到 了 364.5096%。 7Figure 3:空 间 状 态 Exponetial Smoothing交 易 SP500 8require(quantmod) l

28、ibrary(forecast) require(zoo) symbol-c(“GSPC“) symbolData - new.env() getSymbols(symbol, env = symbolData, src = “yahoo“, from=“2007-01-01“, to=“2009-01-01“) sp500 - do.call(cbind,eapply(symbolData, Cl) regime-function(rolling) fit - ets(rolling,damped=FALSE) fcast - forecast(fit, h=1,level=95) if(t

29、ail(rolling,1)fcast$mean)1 else-1 days - 30 r - na.omit(ROC(sp500,1,“discrete“) TS - zoo(r) side - rollapply(TS, width = days, by = 1, FUN = regime, by.column=FALSE, align = “right“) side - as.numeric(head(side,-1) next.r - tail(r,-days) trades - side*next.r layout(1:2) profit - cumprod(trades+1) sp

30、500.plot - tail(sp500,-(days+1) plot(sp500.plot, main = “SP500 Prices“) plot(profit) length(sp500.plot) length(profit) 0.4 百 百 百 度 度 度 文 文 文 库 库 库 资 资 资 源 源 源 1. R语 言 时 间 序 列 中 文教 程 2. R语 言 金 融 工 程 中 文教 程 93. R语 言 样 品 比 较 应 用 举 例 4. R语 言 德 国 坦 克 问 题 5. R语 言 医 药 研 究 6. R语 言 形 成 三 角 形 的 概 率 7. R语 言 小 波 分 析 8. R语 言 太 空 总 署 的 蓄 电 池 10