3.1 基础模型

fGarch : Various submodels arise from this model, and are passed to the ugarchspec “variance.model” list via the submodel option,

• The simple GARCH model of Bollerslev (1986) when \(\lambda = \delta = 2\) and \(\eta_{1j} = \eta_{2j} = 0\) (submodel = ‘GARCH’).
• The Absolute Value GARCH (AVGARCH) model of Taylor (1986) and Schwert (1990) when \(\lambda = \delta = 1\) and \(|\eta_{1j} | ≤ 1\) (submodel = ‘AVGARCH’).
• The GJR GARCH (GJRGARCH) model of Glosten et al. (1993) when \(\lambda = \delta = 2\) and \(\eta_{2j} = 0\) (submodel = ‘GJRGARCH’).
• The Threshold GARCH (TGARCH) model of Zakoian (1994) when \(\lambda = \delta = 1, \eta_{2j} = 0\) and \(|\eta_{1j} | ≤ 1\) (submodel = ‘TGARCH’).
• The Nonlinear ARCH model of Higgins et al. (1992) when \(\delta = \lambda\) and \(\eta_{1j} = \eta_{2j} = 0\) (submodel = ‘NGARCH’).
• The Nonlinear Asymmetric GARCH model of Engle and Ng (1993) when \(\delta = \lambda = 2\) and \(\eta_{1j} = 0\) (submodel = ‘NAGARCH’).
• The Asymmetric Power ARCH model of Ding et al. (1993) when \(\delta = \lambda, \eta_{2j} = 0\) and \(|\eta_{1j} | ≤ 1\) (submodel = ‘APARCH’).
• The Exponential GARCH model of Nelson (1991) when \(\delta = 1, \lambda = 0\) and \(\eta_{2j} = 0\) (not implemented as a submodel of fGARCH).
• The Full fGARCH model of Hentschel (1995) when \(\delta = \lambda\) (submodel = ‘ALLGARCH’).

The choice of distribution is entered via the ‘distribution.model’ option of the ugarchspec method. The package also implements a set of functions to work with the parameters of these distributions. These are:

• ddist(distribution = "norm", y, mu = 0, sigma = 1, lambda = -0.5, skew = 1, shape = 5). The density (d*) function.
• pdist(distribution = "norm", q, mu = 0, sigma = 1, lambda = -0.5, skew = 1, shape = 5). The distribution (p*) function.
• qdist(distribution = "norm", p, mu = 0, sigma = 1, lambda = -0.5, skew = 1, shape = 5). The quantile (q*) function.
• rdist(distribution = "norm", n, mu = 0, sigma = 1, lambda = -0.5, skew = 1, shape = 5). The sampling (q*) function.
• fitdist(distribution = "norm", x, control = list()). A function for fitting data using any of the included distributions.
• dskewness(distribution = "norm", skew = 1, shape = 5, lambda = -0.5). The distribution skewness (analytical where possible else by quadrature integration).
• dkurtosis(distribution = "norm", skew = 1, shape = 5, lambda = -0.5). The distribution excess kurtosis (analytical where it exists else by quadrature integration).

The family of APARCH models includes the ARCH and GARCH models, and five other ARCH extensions as special cases:

• ARCH Model of Engle when \(\delta = 2\), \(\gamma_{i} = 0\), and \(\beta_{j} = 0\).
• GARCH Model of Bollerslev when \(\delta = 2\), and \(\gamma_{i} = 0\).
• TS-GARCH Model of Taylor and Schwert when \(\delta = 1\), and \(\gamma_{i} = 0\).
• GJR-GARCH Model of Glosten, Jagannathan, and Runkle when \(\delta = 2\).
• T-ARCH Model of Zakoian when \(\delta = 1\).
• N-ARCH Model of Higgens and Bera when \(\gamma_{i} = 0\), and \(\beta_{j} = 0\).
• Log-ARCH Model of Geweke and Pentula when \(\delta → 0\).

原文：Parameter Estimation of ARMA Models with GARCH/APARCH Errors - An R and SPlus Software Implementation文献中的2. Mean and Variance Equation。

有关多种GARCH模式比较，请参考参考文献中的链接3… 包括比较：

• auto.arima
• exponential smoothing models (ETS)
• GARCH (包括GARCH、eGARCH、iGARCH、fGARCH、gjrGARCH等模式)
• exponential weighted models

\[\begin{equation} \sigma^2_{t} = \omega + \sum_{i=1}^{\rho}(\alpha_{i} + \gamma_{i} I_{t-i}) \varepsilon_{t-i}^{2} + \sum_{j=1}^{q}\beta_{j}\sigma^{2}_{t-j}\ \cdots\ Equation\ 3.1.1 \end{equation}\]

在之前的文章已经分别比较多种统计模式，得知GJR-GARCH模型的预测结果最为精准，以下稍微介绍下平滑移动加权模型。

3.2 ARMA 模型

ARMA Mean Equation: The ARMA(p,q) process of autoregressive order p and moving average order q can be described as

\[\begin{align*} x_{t} &= \mu + \sum_{i=1}^{m} \alpha_{i} x_{t-i} + \sum^{n}_{j=1} \beta_{j} \varepsilon_{t-j} + \varepsilon_{t} \\ &= \mu + \alpha(B)x_{t} + \beta(B) \varepsilon_{t} \end{align*} \cdots\ Equation\ 3.2.1\]

以上函数乃滑动加权指数，请参阅以下链接以了解更多详情：

• Computer Lab Sessions 2&3
• 時間序列分析 - 總體經濟與財務金融之應用 - 定態時間序列II ARMA模型
• Introduction to the rugarch package
• Parameter Estimation of ARMA Models with GARCH/APARCH Errors - An R and SPlus Software Implementation
• How to choose the order of a GARCH model?

3.3 GJR-GARCH：ARMA(p,q)值最优化（旧程序）

使用极大似然法计算最优arma order中的p值与q值，将原本默认值最优化，让模型马尔可夫化，每日的p与q值最将会使用最佳值… 不包括d值。

armaSearch <- suppressWarnings(function(data, .method = 'CSS-ML') {

    ## I set .method = 'CSS-ML' as default method since the AIC value we got is 
    ##  smaller than using method 'ML' while using method 'CSS' facing error.
    ## 
    ## https://stats.stackexchange.com/questions/209730/fitting-methods-in-arima
    ## According to the documentation, this is how each method fits the model:
    ##  - CSS minimises the sum of squared residuals.
    ##  - ML maximises the log-likelihood function of the ARIMA model.
    ##  - CSS-ML mixes both methods: first, CSS is run, the starting parameters 
    ##    for the optimization algorithm are set to zeros or to the values given 
    ##    in the optional argument init; then, ML is applied passing the CSS 
    ##    parameter estimates as starting parameter values for the optimization algorithm.
    
    .methods = c('CSS-ML', 'ML', 'CSS')
    
    if(!.method %in% .methods) 
      stop(paste('Kindly choose .method among ', 
                 paste0(.methods, collapse = ', '), '!'))
    
    armacoef <- data.frame()
    for (p in 0:5){
      for (q in 0:5) {
        #data.arma = arima(diff(data), order = c(p, 0, q))
        #'@ data.arma = arima(data, order = c(p, 1, q), method = .method)
        if(.method == 'CSS-ML') {
          data.arma = tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'CSS-ML')
            mth = 'CSS-ML'
            list(arma, mth)
          }, error = function(e) tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'ML')
            mth = 'ML'
            list(arma = arma, mth = mth)
          }, error = function(e) {
            arma = arima(data, order = c(p, 1, q), method = 'CSS')
            mth = 'CSS'
            list(arma = arma, mth = mth)
          }))
          
        } else if(.method == 'ML') {
          data.arma = tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'ML')
            mth = 'ML'
            list(arma = arma, mth = mth)
          }, error = function(e) tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'CSS-ML')
            mth = 'CSS-ML'
            list(arma = arma, mth = mth)
          }, error = function(e) {
            arma = arima(data, order = c(p, 1, q), method = 'CSS')
            mth = 'CSS'
            list(arma = arma, mth = mth)
          }))
          
        } else if(.method == 'CSS') {
          data.arma = tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'CSS')
            mth = 'CSS'
            list(arma = arma, mth = mth)
          }, error = function(e) tryCatch({
            arma = arima(data, order = c(p, 1, q), method = 'CSS-ML')
            mth = 'CSS-ML'
            list(arma = arma, mth = mth)
          }, error = function(e) {
            arma = arima(data, order = c(p, 1, q), method = 'ML')
            mth = 'ML'
            list(arma = arma, mth = mth)
          }))
          
        } else {
          stop(paste('Kindly choose .method among ', paste0(.methods, collapse = ', '), '!'))
        }
        names(data.arma) <- c('arma', 'mth')
        
        #cat('p =', p, ', q =', q, 'AIC =', data.arma$arma$aic, '\n')
        armacoef <- rbind(armacoef,c(p, q, data.arma$arma$aic))
      }
    }
    
    ## ARMA Modeling寻找AIC值最小的p,q
    colnames(armacoef) <- c('p', 'q', 'AIC')
    pos <- which(armacoef$AIC == min(armacoef$AIC))
    cat(paste0('method = \'', data.arma$mth, '\', the min AIC = ', armacoef$AIC[pos], 
               ', p = ', armacoef$p[pos], ', q = ', armacoef$q[pos], '\n'))
    return(armacoef)
  })

然后把以上的函数嵌入以下GARCH模型，将原本固定参数的ARMA值浮动化。

calC <- function(mbase, currency = 'JPY=X', ahead = 1, price = 'Cl') {
  
  # Using "memoise" to automatically cache the results
  source('function/filterFX.R')
  source('function/armaSearch.R')
  mbase = suppressWarnings(filterFX(mbase, currency = currency, price = price))
  
  armaOrder = suppressWarnings(armaSearch(mbase))
  armaOrder %<>% dplyr::filter(AIC == min(AIC)) %>% .[c('p', 'q')] %>% unlist
  
  spec = ugarchspec(
    variance.model = list(
      model = 'gjrGARCH', garchOrder = c(1, 1), 
      submodel = NULL, external.regressors = NULL, 
      variance.targeting = FALSE), 
    mean.model = list(
      armaOrder = armaOrder, 
      include.mean = TRUE, archm = FALSE, 
      archpow = 1, arfima = FALSE, 
      ## https://stats.stackexchange.com/questions/73351/how-does-one-specify-arima-p-d-q-in-ugarchspec-for-ugarchfit-in-rugarch?answertab=votes#tab-top
      ## https://d.cosx.org/d/2689-2689/9
      external.regressors = NULL, 
      archex = FALSE), 
    distribution.model = 'snorm')
  fit = ugarchfit(spec, mbase, solver = 'hybrid')
  fc = ugarchforecast(fit, n.ahead = ahead)
  res = tail(attributes(fc)$forecast$seriesFor, 1)
  colnames(res) = names(mbase)
  latestPrice = tail(mbase, 1)

  #rownames(res) <- as.character(forDate)
  latestPrice <- xts(latestPrice)
  #res <- as.xts(res)
  
  tmp = list(latestPrice = latestPrice, forecastPrice = res, 
             AIC = infocriteria(fit))
  return(tmp)
}

3.4 Fi-GJR-GARCH：ARFIMA(p,d,q)值最优化（新程序）

The fractionally integrated GARCH model (‘fiGARCH’) : Contrary to the case of the ARFIMA model, the degree of persistence in the FIGARCH model operates in the oppposite direction, so that as the fractional differencing parameter d gets closer to one, the memory of the FIGARCH process increases, a direct result of the parameter acting on the squared errors rather than the conditional variance. When d = 0 the FIGARCH collapses to the vanilla GARCH model and when d = 1 to the integrated GARCH model…

Motivated by the developments in long memory processes, and in particular the ARFIMA type models (see section 2.1), Baillie et al. (1996) proposed the fractionally integrated generalized autoregressive conditional heteroscedasticity, or FIGARCH, model to capture long memory (in essence hyperbolic memory). Unlike the standard GARCH where shocks decay at an exponential rate, or the integrated GARCH model where shocks persist forever, in the FIGARCH model shocks decay at a slower hyperbolic rate. Consider the standard GARCH equation:

\[\begin{equation} \sigma^{2}_{t} = \omega + \alpha (L) \varepsilon^{2}_{t} + \beta (L) \sigma^{2} \cdots\ Equation\ 3.1.2 \end{equation}\]

原文：Introduction to the rugarch package文献中的2.2.10 The fractionally integrated GARCH model (‘fiGARCH’)

然后计算最优arma order… 也包括d值。

opt_arma <- function(mbase){
  #ARMA Modeling minimum AIC value of `p,d,q`
  fit <- auto.arima(mbase)
  arimaorder(fit)
  }

再来就设置Garch模型中的arfima参数，将原本固定的d值浮动化。

calc_fx <- function(mbase, currency = 'JPY=X', ahead = 1, price = 'Cl') {
  
  ## Using "memoise" to automatically cache the results
  ## http://rpubs.com/englianhu/arma-order-for-garch
  source('function/filterFX.R')
  #'@ source('function/armaSearch.R') #old optimal arma p,q value searching, but no d value. 
  source('function/opt_arma.R') #rename the function best.ARMA()
  
  mbase = suppressWarnings(filterFX(mbase, currency = currency, price = price))
  armaOrder = opt_arma(mbase)
  
  ## Set arma order for `p, d, q` for GARCH model.
  #'@ https://stats.stackexchange.com/questions/73351/how-does-one-specify-arima-p-d-q-in-ugarchspec-for-ugarchfit-in-rugarch
  spec = ugarchspec(
    variance.model = list(
      model = 'gjrGARCH', garchOrder = c(1, 1), 
      submodel = NULL, external.regressors = NULL, 
      variance.targeting = FALSE), 
    mean.model = list(
      armaOrder = armaOrder[c(1, 3)], #set arma order for `p` and `q`.
      include.mean = TRUE, archm = FALSE, 
      archpow = 1, arfima = TRUE, #set arima = TRUE
      external.regressors = NULL, 
      archex = FALSE), 
    fixed.pars = list(arfima = armaOrder[2]), #set fixed.pars for `d` value
    distribution.model = 'snorm')
  
  fit = ugarchfit(spec, mbase, solver = 'hybrid')
  
  fc = ugarchforecast(fit, n.ahead = ahead)
  #res = xts::last(attributes(fc)$forecast$seriesFor)
  res = tail(attributes(fc)$forecast$seriesFor, 1)
  colnames(res) = names(mbase)
  latestPrice = tail(mbase, 1)

  #rownames(res) <- as.character(forDate)
  latestPrice <- xts(latestPrice)
  #res <- as.xts(res)
  
  tmp = list(latestPrice = latestPrice, forecastPrice = res, 
             AIC = infocriteria(fit))
  return(tmp)
  }

4.1 运行时间

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

## 测试运行时间。
#'@ 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]))

##    user  system elapsed 
##   8.240   0.024   8.382

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

## method = 'CSS-ML', the min AIC = -11085.3927792844, p = 5, q = 2

##    user  system elapsed 
##   12.82    0.00   13.61

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

4.2 数据误差率

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

#'@ 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]

#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()函数预测汇价数据。

## ------------- 模拟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()如下。

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')
    }
  }

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

## ------------- 模拟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)
  }

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

## calC()模拟数据误差率
task_progress(.pattern = '^pred1', .loops = FALSE)

## Current Tokyo Time : 2018-09-07 04:27:06 
## 
##      .id    x    n progress
## 1 USDAUD 1214 1215   99.92%
## 2 USDEUR 1214 1215   99.92%
## 3 USDGBP 1215 1216   99.92%
## 4 USDCHF 1214 1215   99.92%
## 5 USDCAD 1214 1214  100.00%
## 6 USDCNY 1214 1215   99.92%
## 7 USDJPY 1213 1215   99.84%
## 
## ================ 99.92% ================

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

## Current Tokyo Time : 2018-09-07 04:27:06 
## 
##      .id    x    n progress
## 1 USDAUD 1215 1215  100.00%
## 2 USDEUR 1215 1215  100.00%
## 3 USDGBP 1216 1216  100.00%
## 4 USDCHF 1215 1215  100.00%
## 5 USDCAD 1214 1214  100.00%
## 6 USDCNY 1212 1215   99.75%
## 7 USDJPY 1215 1215  100.00%
## 
## ================ 99.96% ================

以上结果显示，模拟后的数据的误差率非常渺小¹。以下筛选pred1与pred2同样日期的有效数据。

##数据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)

## 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原始数据有误就不是统计模型的问题了。

fx %<>% dplyr::filter(!Date %in% ntimeID)

4.3 精准度

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

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')

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')

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'))

Group Table Summary

Cat	mse	n	Akaike.mse	Akaike	Bayes	Shibata	Hannan.Quinn
pred1	0.833286	33544	0.0001888	-5.687425	-5.549847	-5.690497	-5.632126
pred2	0.066189	33560	0.0002344	-6.235520	-6.138908	-6.236992	-6.196687

6.1 文件与系统资讯

以下乃此文献资讯：

• 文件建立日期：2018-08-07
• 文件最新更新日期：2018-09-07
• R version 3.4.4 (2018-03-15)
• R语言版本：3.4.4
• rmarkdown 程序包版本：1.10.8
• 文件版本：1.0.1
• 作者简历：®γσ, Eng Lian Hu
• GitHub：源代码
• 其它系统资讯：

System Summary

Category	session_info	Category	Sys.info
version	R version 3.4.4 (2018-03-15)	sysname	Linux
system	x86_64, linux-gnu	release	4.4.0-111-generic
ui	X11	version	#134~14.04.1-Ubuntu SMP Mon Jan 15 15:39:56 UTC 2018
language	(EN)	nodename	4b21bf3ded14
collate	C.UTF-8	machine	x86_64
tz	Etc/UTC	login	unknown
date	2018-09-06	user	rstudio-user
Current time	2018-09-07 04:27:37 JST	effective_user	rstudio-user

6.2 参考文献

1. How does one specify arima (p,d,q) in ugarchspec for ugarchfit in rugarch?
2. How to read p,d and q of auto.arima()?
3. binary.com : Job Application - Quantitative Analyst

---

Powered by - Copyright® Intellectual Property Rights of Scibrokes®個人の経営企業