binary-Q1XREG-ARFIMAX-GARCH.Rmd

---
title: "<img src='www/binary-logo-resize.jpg' width='240'>"
subtitle: "[binary.com](https://猫城.com/englianhu/binary.com-interview-question) 面试试题 I - GARCH模型中的外部加权因素变量参数应用"
author: "[®γσ, Lian Hu](https://englianhu.猫城.io/) <img src='www/RYO.jpg' width='24'> <img src='www/RYU.jpg' width='24'> <img src='www/ENG.jpg' width='24'>®"
date: "`r lubridate::today('Asia/Tokyo')`"
output:
  html_document: 
    number_sections: yes
    toc: yes
    toc_depth: 4
    toc_float:
      collapsed: yes
      smooth_scroll: yes
---

# 简介

在上一篇*GARCH模型中的ARIMA(p,d,q)参数最优化*中将`ARIMA(p,d,q)`值嵌入GARCH模型最优化（至于Garch Order(1,1)则尚未最优化），*binary.com Interview Question I - Interday Betting Strategy Models Comparison (Financial Betting)*采用不同的交易投注模式。今天僕就开始尝试使用`xreg`（外生变量/外部变量/虚拟变量）和`arfimax()`添加外部因素，以形成隐藏马尔可夫模型^[详情请参阅[解密复兴科技 - 基于隐蔽马尔科夫模型的时序分析方法](https://猫城.com/scibrokes/real-time-fxcm/blob/master/reference/%E8%A7%A3%E5%AF%86%E5%A4%8D%E5%85%B4%E7%A7%91%E6%8A%80%20-%20%E5%9F%BA%E4%BA%8E%E9%9A%90%E8%94%BD%E9%A9%AC%E5%B0%94%E7%A7%91%E5%A4%AB%E6%A8%A1%E5%9E%8B%E7%9A%84%E6%97%B6%E5%BA%8F%E5%88%86%E6%9E%90%E6%96%B9%E6%B3%95.pdf)，该书中提及当初复兴科技在隐藏变量中添加宏观因素、估计是国家利率、GDP等因素而导致风险大，后来该公司坚姆斯西门斯大刀阔斧对该模型切割宏观因素。故此得选择不同的外部因素，再经过回测至少3年数据以确保该模型的可行性]。

`rugarch`程序包中的variance.model和mean.model都有`external.regressors`可让僕们自行添加外部因素，由于garch order和arma order的相关多变量函数非常困难，倘若要测试的话得从garch.order(1,1)测试到garch.order(5,5)^[就如[binary.com 面试试题 I - GARCH模型中的ARIMA(p,d,q)参数最优化](http://rpubs.com/englianhu/binary-Q1FiGJRGARCH)中的`armaSearch()`函数。]，欲知更多详情，请参阅[binary.com Interview Question I - Comparison of Univariate GARCH Models](http://rpubs.com/englianhu/binary-Q1Uni-GARCH)中的**7 GARCH(1,1)?**。

```{r setup}
suppressPackageStartupMessages(require('BBmisc'))

## 读取程序包
pkg <- c('lubridate', 'plyr', 'dplyr', 'magrittr', 'stringr', 'rugarch', 'forecast', 'quantmod', 'microbenchmark', 'knitr', 'kableExtra', 'formattable', 'quantmod', 'TTR', 'memoise')
suppressAll(lib(pkg))
rm(pkg)

funs <- c('calc_fx.R', 'opt_arma.R', 'filterFX.R')
l_ply(funs, function(x) source(paste0('./function/', x)))
```

# 数据

首先读取[Binary.com Interview Q1 (Extention)](http://rpubs.com/englianhu/binary-Q1E)的汇市数据。

```{r read-data, warning=FALSE}
cr_code <- c('AUDUSD=X', 'EURUSD=X', 'GBPUSD=X', 'CHF=X', 'CAD=X', 
             'CNY=X', 'JPY=X')

#'@ names(cr_code) <- c('AUDUSD', 'EURUSD', 'GBPUSD', 'USDCHF', 'USDCAD', 
#'@                     'USDCNY', 'USDJPY')

names(cr_code) <- c('USDAUD', 'USDEUR', 'USDGBP', 'USDCHF', 'USDCAD', 'USDCNY', 'USDJPY')

price_type <- c('Op', 'Hi', 'Lo', 'Cl')

## 读取雅虎数据。
mbase <- sapply(names(cr_code), function(x) readRDS(paste0('/home/englianhu/Documents/猫城/binary.com-interview-question-data/data/', x, '.rds')) %>% na.omit)
```

数据简介报告。

```{r data-summary}
sapply(mbase, summary) %>% 
  kable %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>%
  scroll_box(width = '100%', height = '400px')
```

*桌面2.1：数据简介。*

# 统计建模

## 基础gjrGARCH模型

$$
\begin{equation}
σ^2_{t} = \left ( \omega + \sum^{m}_{j=1} \zeta_{j} \nu_{jt} \right ) + \sum^{q}_{j=1} \left ( \alpha_{j}\varepsilon ^{2}_{t-1} + \gamma_{j}I_{t-j}\varepsilon ^{2}_{t-1}\right ) + \sum^{p}_{j=1} \beta_{j}\sigma^{2}_{t-j}
\ \cdots\ Equation\ 3.1.1
\end{equation}
$$

请参阅*Introduction to the `rugarch` Package*以了解更多详情。

## Fi-GJR-GARCH(p,d,q)值最优化

请参阅*GARCH模型中的ARIMA(p,d,q)参数最优化*。

[R学习笔记（2）——SARIMA模型](https://blog.csdn.net/qqqqq1993qqqqq/article/details/51818718)和[`sarima`程序包](https://cran.r-project.org/web/packages/sarima/sarima.pdf)或许能基于季节性因素把原本的ARIMA(p,d,q)值，通过SARIMA(p,d,q)值更为精准。

## FIGARCH

## HYGARCH

# 加权因子

## 外部因子

[`rugarch`包与R语言中的garch族模型](http://www.itread01.com/content/1521893044.html)中提及外生变量。 [`rugarch`:GARCH external regressors](https://quant.stackexchange.com/questions/27491/rugarch-garch-external-regressors)和[Rugarch package using external regressors](https://stackoverflow.com/questions/39193250/rugarch-package-using-external-regressors)和[Different external regressor in rugarch give the same result](http://r.789695.n4.nabble.com/Different-external-regressor-in-rugarch-give-the-same-result-td4687297.html)三篇文章中也有提及`setbounds(spec)<-list(vxreg1=c(-1,1))`。

>And DO NOT USE c(model="fGARCH", submodel="GARCH"). Use directly model="sGARCH".

*原文 : Different external regressor in rugarch give the same result*

不晓得该作者建议使用非fGARCH的目的为何？不过[binary.com Interview Question I](http://rpubs.com/englianhu/binary-Q1)和[binary.com Interview Question I - Comparison of Univariate GARCH Models](http://rpubs.com/englianhu/binary-Q1Uni-GARCH)比较多个GARCH模型，而使用fGARCH添加子模型频频出现不知名的代码错误、而且未使用fGARCH会比较精准。

[rugarch包编程程序](http://bbs.pinggu.org/thread-5533864-1-1.html)也应用xreg，而外生变量是`时间`和`国家`变量。

[R软件如何在GARCH模型中加入虚拟变量？](http://bbs.pinggu.org/thread-2131656-1-1.html)中也有使用虚拟变量（Dummy Variable）。

[rgarch预测股票汇报的疑问](https://d.cosx.org/d/107102-107102)将股票红利变量因素植入GARCH模型中。

```
Hi all, I'm using the rugarch package (which is btw an excellent piece of software). I'm using the A-PARCH(1,1) model with an external regressor (the Yang-Zhang volatility proxy) and I've a doubt about how to perform forecasting using the model: 

  * Suppose I've n observations for the returns of a financial serie and for the external regressor
  * Then, I use ugarchspec with an A-PARCH(0,1) model with the external regressor, the ugarchspec and ugarchfit steps work pretty well.
  * Now I want to do forecasting. I don't quite fully understand the limitation of this process. When I use ugarchforecast, the main limitation is that I only can forecast 1-step ahead due to the external regressor data limitation (I don't have more than n data for the external regressor)??? 

Now, if I fit a model with the n external regressor data and if I do a 1-step ahead forecast for t+1, then the next day I recieve new data for the external regressor, how to use the rugarch package to use this new data and generate a new 1-step ahead forecast from t+1 to t+2? is this possible? 

Thanks in advance, I've read a lot of the rugarch documentation and from this mailing list, but I'm a bit lost regarding to my doubts. 

Any help will be really appreciated. 

Best, 


-- 
Diego Ignacio Acuña Rozas 
Mg. (c) Ciencias de la Ingeniería Informática 
Universidad Técnica Federico Santa María 
#==================================================================
As this appears to be a general point of confusion, see below a commented example: 
######################### 
library(rugarch) 
library(xts) 
data(sp500ret) 
spx<-xts(sp500ret, as.Date(rownames(sp500ret))) 
xreg<-xts(rnorm(nrow(spx)), index(spx)) 
colnames(xreg)<-"xreg" 
# assume xreg is an actual series, so we lag it 
# as we would do in a real application 
xreg = lag(xreg,1) 
inputs<-na.omit(cbind(spx, xreg, join="left")) 
# real time forecasting 
spec<-ugarchspec(mean.model=list(external.regressors=inputs[1:2000,2])) 
fit<-ugarchfit(spec, inputs[1:2000,1]) 
# 2 ways to do real-time forecasting (ugarchforecast and ugarchfilter) 
# Example: forecast points 2001:2020 
xforc = xts(matrix(NA, ncol=2, nrow=20), index(inputs[2001:2020])) 
sforc = xts(matrix(NA, ncol=2, nrow=20), index(inputs[2001:2020])) 
for(i in 1:20){ 
   # Forecast(T+1)|Information(T) 
   # 1. Create a similar spec as you used in estimation 
   # and add the lagged regressor upto time T 
specf1<-ugarchspec(mean.model=list(external.regressors=inputs[1:(2000+i-1),2])) 
   # Pass the estimated coefficients from the estimation upto time 2000 
   setfixed(specf1)<-as.list(coef(fit)) 
   # 2. Forecast using ugarchforecast on a specification with fixed 
parameters 
   # where n.old is used in order to recreate the correct start-up 
conditions 
   # used in the fitting routine 
   f1<-ugarchforecast(specf1, inputs[1:(2000+i-1),1], n.ahead=1, n.old=2000) 
   # 3. Forecast using ugarchfilter on a specification with fixed 
parameters. 
   # For this method, append a new row to the end of the data with zeros, 
   # as you would do with related filters. This forces the routine to 
   # output the value at time T+1 
   newdat<-rbind(inputs[1:(2000+i-1),],xts(matrix(0, nrow=1, ncol=2), 
tail(move(index(inputs[1:(2000+i-1)])),1))) 
   specf2<-ugarchspec(mean.model=list(external.regressors=newdat[,2])) 
   setfixed(specf2)<-as.list(coef(fit)) 
   f2<-ugarchfilter(specf2, newdat[,1], n.old=2000) 
   # fitted = estimated conditional mean values for uGARCHfit objects 
   # fitted = forecast/filtered conditional mean values for 
uGARCHforecast/uGARCHfilter objects 
   xforc[i,1] = as.numeric(fitted(f1)) 
   xforc[i,2] = as.numeric(tail(fitted(f2),1)) 
   # sigma = conditional sigma values (fitted/forecast etc) 
   sforc[i,1] = as.numeric(sigma(f1)) 
   sforc[i,2] = as.numeric(tail(sigma(f2),1)) 
} 
# check 
all.equal(xforc[,1], xforc[,2]) 
all.equal(sforc[,1], sforc[,2]) 
# check that the 1-ahead forecast directly from the fitted object is also 
# the same 
all.equal(as.numeric(xforc[1,1]), as.numeric(fitted(ugarchforecast(fit, 
n.ahead=1)))) 
all.equal(as.numeric(sforc[1,1]), as.numeric(sigma(ugarchforecast(fit, 
n.ahead=1)))) 
# check the filter values vs the fitted values (i.e. why we use the 
n.old argument) 
all.equal(fitted(fit), fitted(f2)[1:2000]) 
all.equal(sigma(fit), sigma(f2)[1:2000]) 
######################### 


Regards, 

Alexios 
```

以上的文章，[Rugarch package using external regressors, a forecasting doubt](http://r.789695.n4.nabble.com/Rugarch-package-using-external-regressors-a-forecasting-doubt-td4709073.html)举例，使用`rnorm()`加权因子和。

> Note that there does not seem to be an option to use SARMA models in the "rugarch" package, so you will have to let the "S" part go. But if there is a seasonal pattern (and that is quite likely when it comes to tourist arrivals), you will have to account for it somehow. Consider using exogenous seasonal variables (dummies or Fourier terms) in the conditional mean model via the argument external.regressors inside the argument mean.model in function ugarchspec. Alternatively, note that a SARMA model corresponds to a restricted ARMA model. An approximation of SARMA could thus be an ARMA with the appropriate lag order but without the SARMA-specific parameter restrictions (since those might not be available in "rugarch").

以上文章提及季节因素，[Fitting ARIMA-GARCH model using “rugarch” package](https://stats.stackexchange.com/questions/176550/fitting-arima-garch-model-using-rugarch-package?answertab=votes#tab-top)。欲建模`sarima`可使用[`sarima`程序包](https://cran.r-project.org/web/packages/sarima/sarima.pdf)

请参阅*Introduction to the `rugarch` Package*以了解更多详情。

## ROC函数

我们可以使用`ROC()`计算$t$和$t-1$的价格是否高于

[ROC和AUC介绍以及如何计算AUC](http://alexkong.net/2013/06/introduction-to-auc-and-roc/)介绍ROC（Receiver Operating Characteristic）曲线和AUC（Area Under Curve）常被用来评价一个二值分类器（binary classifier）的优劣。该篇博文简单介绍ROC和AUC的特点，以及更为深入地，讨论如何作出ROC曲线图以及计算AUC。

[《金融时间序列预测》：第十二章：鄀量化投资初步](http://blog.sciencenet.cn/blog-577790-830021.html)就`quantmod`包的回测方法介绍基于技术指标的量化投资方法。

[第一次使用鄀語言做回測：六分鐘，就上手!](http://www.bituzi.com/2014/12/Rbacktest6mins.html)

[[原]量化投资教程：用R语言打造量化分析Web平台](https://segmentfault.com/a/1190000004543727)介绍了许多交易实用的函数。

[An Example Of A Trading Strategy Coded Using Quantmod Package In R](https://www.quantinsti.com/blog/an-example-of-a-trading-strategy-coded-in-r/)

```{r roc}

```

## ARFIMAX模型



请参阅*Introduction to the `rugarch` Package*以了解更多详情。

# 模式比较

## 运行时间

首先比较运行时间，哪个比较高效。

```{r processing-time, warning = FALSE}
## 测试运行时间。
#'@ microbenchmark(fit <- calc_fx(mbase[[names(cr_code)[sp]]], currency = cr_code[sp]))
#'@ microbenchmark(fit2 <- calC(mbase[[names(cr_code)[sp]]], currency = cr_code[sp]))

## 随机抽样货币数据，测试运行时间。
sp <- sample(1:7, 1)

system.time(fit1 <- calc_fx(mbase[[names(cr_code)[sp]]], currency = cr_code[sp]))
system.time(fit2 <- calC(mbase[[names(cr_code)[sp]]], currency = cr_code[sp]))
```

由于使用`microbenchmark`非常耗时，而且双方实力悬殊，故此僕使用`system.time()`比较运行速度，结果还是新程序`calc_fx()`比旧程序`calC()`迅速。

## 数据误差率

以下僕运行数据测试后事先储存，然后直接读取。首先过滤`timeID`时间参数，然后才模拟预测汇价。

```{r tidy-data}
#'@ ldply(mbase, function(x) range(index(x)))
#     .id         V1         V2
#1 USDAUD 2012-01-02 2017-08-30
#2 USDEUR 2012-01-02 2017-08-30
#3 USDGBP 2012-01-02 2017-08-30
#4 USDCHF 2012-01-02 2017-08-30
#5 USDCAD 2012-01-02 2017-08-30
#6 USDCNY 2012-01-02 2017-08-30
#7 USDJPY 2012-01-02 2017-08-30

timeID <- llply(mbase, function(x) as.character(index(x))) %>% 
  unlist %>% unique %>% as.Date %>% sort
timeID <- c(timeID, xts::last(timeID) + days(1)) #the last date + 1 in order to predict the next day of last date to make whole dataset completed.
timeID0 <- ymd('2013-01-01')
timeID <- timeID[timeID >= timeID0]

## ---------------- 6个R进程并行运作 --------------------
start <- seq(1, length(timeID), ceiling(length(timeID)/6))
#[1]    1  204  407  610  813 1016

stop <- c((start - 1)[-1], length(timeID))
#[1]  203  406  609  812 1015 1217

cat(paste0('\ntimeID <- timeID[', paste0(start, ':', stop), ']'), '\n')
#timeID <- timeID[1:203]
#timeID <- timeID[204:406]
#timeID <- timeID[407:609]
#timeID <- timeID[610:812]
#timeID <- timeID[813:1015]
#timeID <- timeID[1016:1217]

## Some currency data doesn't open market in speficic date.
#Error:
#data/fx/USDCNY/pred1.2015-04-15.rds saved! #only USDJPY need to review
#data/fx/USDGBP/pred1.2015-12-07.rds saved! #only USDCHF need to review
#data/fx/USDCAD/pred1.2016-08-30.rds saved! #only USDCNY need to review
#data/fx/USDAUD/pred1.2016-11-30.rds saved! #only USDEUR need to review
#data/fx/USDCNY/pred1.2017-01-12.rds saved! #only USDJPY need to review
#data/fx/USDEUR/pred1.2017-02-09.rds saved! #only USDGBP need to review
#timeID <- timeID[timeID > ymd('2017-03-08')]

#data/fx/USDCAD/pred2.2015-06-09.rds saved! #only USDCNY need to review
#data/fx/USDCAD/pred2.2015-06-16.rds saved! #only USDCNY need to review
#data/fx/USDCAD/pred2.2015-06-17.rds saved! #only USDCNY need to review
```

模拟`calC()`函数预测汇价数据。

```{r sim-pred1, eval = FALSE, warning = FALSE}
## ------------- 模拟calC()预测汇价 ----------------------
pred1 <- list()

for (dt in timeID) {
  
  for (i in seq(cr_code)) {
    
    smp <- mbase[[names(cr_code)[i]]]
    dtr <- xts::last(index(smp[index(smp) < dt]), 1) #tail(..., 1)
    smp <- smp[paste0(dtr %m-% years(1), '/', dtr)]
    
    pred1[[i]] <- ldply(price_type, function(y) {
      df = calC(smp, currency = cr_code[i], price = y)
      df = data.frame(Date = index(df[[1]][1]), 
                      Type = paste0(names(df[[1]]), '.', y), 
                      df[[1]], df[[2]], t(df[[3]]))
      names(df)[4] %<>% str_replace_all('1', 'T+1')
      df
    })
    
    if (!dir.exists(paste0('data/fx/', names(pred1[[i]])[3]))) 
      dir.create(paste0('data/fx/', names(pred1[[i]])[3]))
    
    saveRDS(pred1[[i]], paste0(
      'data/fx/', names(pred1[[i]])[3], '/pred1.', 
      unique(pred1[[i]]$Date), '.rds'))
    
    cat(paste0(
      'data/fx/', names(pred1[[i]])[3], '/pred1.', 
      unique(pred1[[i]]$Date), '.rds saved!\n'))
    
    }; rm(i)
  }
```

查询模拟测试进度的函数`task_progress()`如下。

```{r check-progress}
task_progress <- function(scs = 60, .pattern = '^pred1', .loops = TRUE) {
  ## ------------- 定时查询进度 ----------------------
  ## 每分钟自动查询与更新以上模拟calC()预测汇价进度（储存文件量）。
  
  if (.loops == TRUE) {
    while(1) {
      cat('Current Tokyo Time :', as.character(now('Asia/Tokyo')), '\n\n')
      
      z <- ldply(mbase, function(dtm) {
        y = index(dtm)
        y = y[y >= timeID0]
        
        cr = as.character(unique(substr(names(dtm), 1, 6)))
        x = list.files(paste0('./data/fx/', cr), pattern = .pattern) %>% 
          str_extract_all('[0-9]{4}-[0-9]{2}-[0-9]{2}') %>% 
          unlist %>% as.Date %>% sort
        x = x[x >= y[1] & x <= xts::last(y)]
        
        data.frame(.id = cr, x = length(x), n = length(y)) %>% 
        mutate(progress = percent(x/n))
      })# %>% tbl_df
      
      print(z)
      
      prg = sum(z$x)/sum(z$n)
      cat('\n================', as.character(percent(prg)), '================\n\n')
      
      if (prg == 1) break #倘若进度达到100%就停止更新。
      
      Sys.sleep(scs) #以上ldply()耗时3~5秒，而休息时间60秒。
    }
  } else {
    
    cat('Current Tokyo Time :', as.character(now('Asia/Tokyo')), '\n\n')
      
    z <- ldply(mbase, function(dtm) {
      y = index(dtm)
      y = y[y >= timeID0]
      
      cr = as.character(unique(substr(names(dtm), 1, 6)))
      x = list.files(paste0('./data/fx/', cr), pattern = .pattern) %>% 
          str_extract_all('[0-9]{4}-[0-9]{2}-[0-9]{2}') %>% 
          unlist %>% as.Date %>% sort
      x = x[x >= y[1] & x <= xts::last(y)]
      
      data.frame(.id = cr, x = length(x), n = length(y)) %>% 
        mutate(progress = percent(x/n))
      })# %>% tbl_df
    
    print(z)
    
    prg = sum(z$x)/sum(z$n)
    cat('\n================', as.character(percent(prg)), '================\n\n')
    }
  }
```

```{r check-files, echo = FALSE, eval = FALSE}
## ------------- 查询缺失文件 ----------------------
## 查询缺失文件。
dts <- sapply(mbase, function(x) {
  y = index(x)
  y[y >= timeID0]
  })

sapply(mbase, function(x) as.character(index(x)) %>% as.Date %>% sort)

fls <- sapply(names(cr_code), function(x) {
   list.files(paste0('./data/fx/', x), pattern = '^pred1') %>% 
     str_extract_all('[0-9]{4}-[0-9]{2}-[0-9]{2}') %>% 
	 unlist %>% as.Date %>% sort
   })

sapply(fls, function(x) timeID[!timeID %in% x] %>% sort)

timeID <- llply(fls, function(x) timeID[!timeID %in% x] %>% sort) %>% unlist %>% as.Date %>% sort
names(timeID) <- NULL
timeID %<>% unique
```

模拟`calc_fx()`函数预测汇价数据。

```{r sim-pred2, eval = FALSE, warning = FALSE}
## ------------- 模拟calc_fx()预测汇价 ----------------------
pred2 <- list()

for (dt in timeID) {
  
  for (i in seq(cr_code)) {
    
    smp <- mbase[[names(cr_code)[i]]]
    dtr <- xts::last(index(smp[index(smp) < dt]), 1) #tail(..., 1)
    smp <- smp[paste0(dtr %m-% years(1), '/', dtr)]
    
    pred2[[i]] <- ldply(price_type, function(y) {
      df = calc_fx(smp, currency = cr_code[i], price = y)
      df = data.frame(Date = index(df[[1]][1]), 
                      Type = paste0(names(df[[1]]), '.', y), 
                      df[[1]], df[[2]], t(df[[3]]))
      names(df)[4] %<>% str_replace_all('1', 'T+1')
      df
    })
    
    if (!dir.exists(paste0('data/fx/', names(pred2[[i]])[3]))) 
      dir.create(paste0('data/fx/', names(pred2[[i]])[3]))
    
    saveRDS(pred2[[i]], paste0(
      'data/fx/', names(pred2[[i]])[3], '/pred2.', 
      unique(pred2[[i]]$Date), '.rds'))
    
    cat(paste0(
      'data/fx/', names(pred2[[i]])[3], '/pred2.', 
      unique(pred2[[i]]$Date), '.rds saved!\n'))
    
    }; rm(i)
  }
```

模拟完毕后，再来就查看数据结果。

```{r data-error}
## calC()模拟数据误差率
task_progress(.pattern = '^pred1', .loops = FALSE)

## calc_fx()模拟数据误差率
task_progress(.pattern = '^pred2', .loops = FALSE)
```

以上结果显示，模拟后的数据的误差率非常渺小^[一些数据模拟时，出现不知名错误。]。以下筛选`pred1`与`pred2`同样日期的有效数据。

```{r tidy-data2}
##数据1
fx1 <- llply(names(cr_code), function(x) {
    fls <- list.files(paste0('data/fx/', x), pattern = '^pred1')
    dfm <- ldply(fls, function(y) {
        readRDS(paste0('data/fx/', x, '/', y))
    }) %>% data.frame(Cat = 'pred1', .) %>% tbl_df
    names(dfm)[4:5] <- c('Price', 'Price.T1')
    dfm
 })
names(fx1) <- names(cr_code)

##数据2
fx2 <- llply(names(cr_code), function(x) {
    fls <- list.files(paste0('data/fx/', x), pattern = '^pred2')
    dfm <- ldply(fls, function(y) {
        readRDS(paste0('data/fx/', x, '/', y))
    }) %>% data.frame(Cat = 'pred2', .) %>% tbl_df
    names(dfm)[4:5] <- c('Price', 'Price.T1')
    dfm
 })
names(fx2) <- names(cr_code)

#合并，并且整理数据。
fx1 %<>% ldply %>% tbl_df
fx2 %<>% ldply %>% tbl_df
fx <- suppressAll(bind_rows(fx1, fx2) %>% arrange(Date) %>% 
  mutate(.id = factor(.id), Cat = factor(Cat)) %>% 
  ddply(.(Cat, Type), function(x) {
    x %>% mutate(Price.T1 = lag(Price.T1, 1))
  }) %>% tbl_df %>% 
    dplyr::filter(Date >= ymd('2013-01-01') & Date <= ymd('2017-08-30')))

rm(fx1, fx2)
```

```{r tidy-data3}
## filter all predictive error where sd >= 20%.
notID <- fx %>% mutate(diff = abs(Price.T1/Price), se = ifelse(diff <= 0.8 | diff >= 1.25, 1, 0)) %>% dplyr::filter(se == 1)
ntimeID <- notID %>% .$Date %>% unique
notID %>% 
  kable(caption = 'Error data') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>%
  scroll_box(width = '100%', height = '400px')
```

僕尝试运行好几次，`USDCHF`都是获得同样的结果。然后将默认的`snorm`分布更换为`norm`就没有出现错误。至于`USDCNY`原始数据有误就不是统计模型的问题了。

```{r tidy-data4}
fx %<>% dplyr::filter(!Date %in% ntimeID)
```

## 精准度

现在就比较下双方的MSE值与AIC值。

```{r aic1}
acc <- ddply(fx, .(Cat, Type), summarise, 
             mse = mean((Price.T1 - Price)^2), 
             n = length(Price), 
             Akaike.mse = (-2*mse)/n+2*4/n, 
             Akaike = mean(Akaike), 
             Bayes = mean(Bayes), 
             Shibata = mean(Shibata), 
             Hannan.Quinn = mean(Hannan.Quinn)) %>% 
  tbl_df %>% mutate(mse = round(mse, 6)) %>% 
  arrange(Type)

acc %>% 
  kable(caption = 'Group Table Summary') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>% 
  group_rows('USD/AUD Open', 1, 2, label_row_css = 'background-color: #e68a00; color: #fff;') %>%
  group_rows('USD/AUD High', 3, 4, label_row_css = 'background-color: #e68a00; color: #fff;') %>%
  group_rows('USD/AUD Low', 5, 6, label_row_css = 'background-color: #e68a00; color: #fff;') %>%
  group_rows('USD/AUD Close', 7, 8, label_row_css = 'background-color: #e68a00; color: #fff;') %>%
  group_rows('USD/EUR Open', 9, 10, label_row_css = 'background-color: #6666ff; color: #fff;') %>%
  group_rows('USD/EUR High', 11, 12, label_row_css = 'background-color: #6666ff; color: #fff;') %>%
  group_rows('USD/EUR Low', 13, 14, label_row_css = 'background-color:#6666ff; color: #fff;') %>%
  group_rows('USD/EUR Close', 15, 16, label_row_css = 'background-color: #6666ff; color: #fff;') %>%
  group_rows('USD/GBP Open', 17, 18, label_row_css = 'background-color: #339966; color: #fff;') %>%
  group_rows('USD/GBP High', 19, 20, label_row_css = 'background-color: #339966; color: #fff;') %>%
  group_rows('USD/GBP Low', 21, 22, label_row_css = 'background-color: #339966; color: #fff;') %>%
  group_rows('USD/GBP Close', 23, 24, label_row_css = 'background-color: #339966; color: #fff;') %>%
  group_rows('USD/CHF Open', 25, 26, label_row_css = 'background-color: #808000; color: #fff;') %>%
  group_rows('USD/CHF High', 27, 28, label_row_css = 'background-color: #808000; color: #fff;') %>%
  group_rows('USD/CHF Low', 29, 30, label_row_css = 'background-color: #808000; color: #fff;') %>%
  group_rows('USD/CHF Close', 31, 32, label_row_css = 'background-color: #808000; color: #fff;') %>%
  group_rows('USD/CAD Open', 33, 34, label_row_css = 'background-color: #666; color: #fff;') %>%
  group_rows('USD/CAD High', 35, 36, label_row_css = 'background-color: #666; color: #fff;') %>%
  group_rows('USD/CAD Low', 37, 38, label_row_css = 'background-color: #666; color: #fff;') %>%
  group_rows('USD/CAD Close', 39, 40, label_row_css = 'background-color: #666; color: #fff;') %>%
  group_rows('USD/CNY Open', 41, 42, label_row_css = 'background-color: #e60000; color: #fff;') %>%
  group_rows('USD/CNY High', 43, 44, label_row_css = 'background-color: #e60000; color: #fff;') %>%
  group_rows('USD/CNY Low', 45, 46, label_row_css = 'background-color: #e60000; color: #fff;') %>%
  group_rows('USD/CNY Close', 47, 48, label_row_css = 'background-color: #e60000; color: #fff;') %>%
  group_rows('USD/JPY Open', 49, 50, label_row_css = 'background-color: #ff3377; color: #fff;') %>%
  group_rows('USD/JPY High', 51, 52, label_row_css = 'background-color: #ff3377; color: #fff;') %>%
  group_rows('USD/JPY Low', 53, 54, label_row_css = 'background-color: #ff3377; color: #fff;') %>%
  group_rows('USD/JPY Close', 55, 56, label_row_css = 'background-color: #ff3377; color: #fff;') %>%
  scroll_box(width = '100%', height = '400px')
```

```{r aic2}
acc <- ddply(fx, .(Cat, .id), summarise, 
             mse = mean((Price.T1 - Price)^2), 
             n = length(Price), 
             Akaike.mse = (-2*mse)/n+2*4/n, 
             Akaike = mean(Akaike), 
             Bayes = mean(Bayes), 
             Shibata = mean(Shibata), 
             Hannan.Quinn = mean(Hannan.Quinn)) %>% 
  tbl_df %>% mutate(mse = round(mse, 6)) %>% 
  arrange(.id)

acc %>% 
  kable(caption = 'Group Table Summary') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive')) %>%
  group_rows('USD/AUD', 1, 2, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/CAD', 3, 4, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/CHF', 5, 6, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/CNY', 7, 8, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/EUR', 9, 10, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/GBP', 11, 12, label_row_css = 'background-color: #003399; color: #fff;') %>%
  group_rows('USD/JPY', 13, 14, label_row_css = 'background-color: #003399; color: #fff;') %>% 
  scroll_box(width = '100%', height = '400px')
```

```{r aic3}
acc <- ddply(fx, .(Cat), summarise, 
             mse = mean((Price.T1 - Price)^2), 
             n = length(Price), 
             Akaike.mse = (-2*mse)/n+2*4/n, 
             Akaike = mean(Akaike), 
             Bayes = mean(Bayes), 
             Shibata = mean(Shibata), 
             Hannan.Quinn = mean(Hannan.Quinn)) %>% 
  tbl_df %>% mutate(mse = round(mse, 6))

acc %>% 
  kable(caption = 'Group Table Summary') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive'))
```

# 结论

# 附录

## 文件与系统资讯

以下乃此文献资讯：

- 文件建立日期：2018-09-07
- 文件最新更新日期：`r today('Asia/Tokyo')`
- `r R.version.string`
- R语言版本：`r getRversion()`
- [**rmarkdown** 程序包](https://猫城.com/rstudio/rmarkdown)版本：`r packageVersion('rmarkdown')`
- 文件版本：1.0.1
- 作者简历：[®γσ, Eng Lian Hu](https://beta.rstudioconnect.com/content/3091/ryo-eng.html)
- 猫城：[源代码](https://猫城.com/englianhu/binary.com-interview-question)
- 其它系统资讯：

```{r info, echo = FALSE, warning = FALSE, results = 'asis'}
suppressMessages(require('dplyr', quietly = TRUE))
suppressMessages(require('formattable', quietly = TRUE))
suppressMessages(require('knitr', quietly = TRUE))
suppressMessages(require('kableExtra', quietly = TRUE))

sys1 <- devtools::session_info()$platform %>% 
  unlist %>% data.frame(Category = names(.), session_info = .)
rownames(sys1) <- NULL

sys1 %<>% rbind(., data.frame(
  Category = 'Current time', 
  session_info = paste(as.character(lubridate::now('Asia/Tokyo')), 'JST'))) %>% 
  dplyr::filter(Category != 'os')

sys2 <- data.frame(Sys.info()) %>% mutate(Category = rownames(.)) %>% .[2:1]
names(sys2)[2] <- c('Sys.info')
rownames(sys2) <- NULL

cbind(sys1, sys2) %>% 
  kable(caption = 'Additional session information:') %>% 
  kable_styling(bootstrap_options = c('striped', 'hover', 'condensed', 'responsive'))

rm(sys1, sys2)
```

## 参考文献

01. [The `rmgarch` Models - Background and Properties](https://raw.猫城usercontent.com/englianhu/binary.com-interview-question/master/reference/The%20rmgarch%20Models%20-%20Background%20and%20Properties.pdf)<img src='www/hot.jpg' width='20'>
02. [binary.com : Job Application - Quantitative Analyst](https://猫城.com/englianhu/binary.com-interview-question)
03. [Introduction to the `rugarch` Package](https://猫城.com/englianhu/binary.com-interview-question/blob/master/reference/Introduction%20to%20the%20rugarch%20Package.pdf)<img src='www/hot.jpg' width='20'>
04. [GARCH模型中的ARIMA(p,d,q)参数最优化](http://rpubs.com/englianhu/binary-Q1FiGJRGARCH)
05. [binary.com Interview Question I - Interday Betting Strategy Models Comparison (Financial Betting)](http://rpubs.com/englianhu/binary-Q1BET)
06. [R语言中ARMA，ARIMA（BOX-JENKINS），SARIMA和ARIMAX模型用于预测时间序列数据](http://tecdat.cn/r%e8%af%ad%e8%a8%80%e4%b8%adarma%ef%bc%8carima%ef%bc%88box-jenkins%ef%bc%89%ef%bc%8csarima%e5%92%8carimax%e6%a8%a1%e5%9e%8b%e7%94%a8%e4%ba%8e%e9%a2%84%e6%b5%8b%e6%97%b6%e9%97%b4%e5%ba%8f%e5%88%97)

---

[<img src="诸子百家考工记/世博量化.png" height="14"/> Sςιβrοκεrs Trαdιηg®](http://www.scibrokes.com)<br>
[**[<img src="诸子百家考工记/世博量化.png" height="14"/> 世博量化®](http://www.scibrokes.com)企业知识产权及版权所有，盗版必究。**]{style="color:RoyalBlue"}