最近接触到了带傅里叶修正项的ARIMA模型（ARIMA with fourier modification, 以下简称FARIMA），学校里面虽然上过一些时间序列，知道time series decomposition 和 ARIMA，但是对seasonal ARIMA以及FARIMA并不是很了解。于是搜集各方资料摸索了一遍用R语言结合FARIMA来进行预测的流程，在这里整理贴出~

1 原理

1.1 ARIMA model

全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA)，是由博克思(Box)和詹金斯(Jenkins)于70年代初提出一著名时间序列预测方法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归， p为自回归项； MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。所谓ARIMA模型，是指将非平稳时间序列转化为平稳时间序列，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。
关于seasonality: 一个时间序列可能有季节性波动的特点，也可能没有。时间序列中的季节性是指数据在S个时期中重复规律变化的模式。对于有季节性的时间序列，在一个周期中其期望值往往是不恒定的，并且会满足以下规律：

因此，在大多数情况下，季节性时间序列是非平稳的。

1.2 Non-seasonal ARIMA model

非非季节性ARIMA模型通常具有ARIMA（p，d，q）形式，其中：
•p是在预测方程中，变量滞后项的数目，称为自回归参数；
•d 是使时间序列成为平稳时所做的差分次数；
•q 是预测方程中残差的滞后项数目，称为移动平均参数。

1.3 Seasonal ARIMA model

对于季节性时间序列，可以通过针对季节性差分使序列平稳，季节性差分是指针对序列值和滞后它S的倍数的序列值的差分。季节性ARIMA模型在multiplicative model中包含非季节性和季节性因子。 SARIMA（p，d，q）（P，D，Q）S的形式，其中：
•p，d，q是如上所述的非季节性ARIMA模型中的参数。
•P是季节性自回归的阶数；
•D是季节性差分的次数；
•Q是季节性移动平均线的阶数；
•S是重复季节性模式的时间跨度。

针对Seasonal ARIMA，我们将数据中的相关性分为两种：regular dependence，指普通的序列间隔导致的相关性；以及seasonal dependence，指间隔S期的观察值之间的相关性。

将季节性囊括进ARIMA有两种方式：
1. 还是将Seasonal ARIMA model按照ARIMA的表达式形式写出，在AR或MA项中加入s阶的滞后算子，但这样会导致项数大大增加。e.g.对于 S=12的月度数据，
如果数据和三年前相同月份的数据相关，那么AR 或MA的阶数就会增加到36。（这里不是很懂，有热心的朋友看懂了希望能帮我解释一下嘻嘻）
2.更常见、简单的一种方式是，将regular dependence与seasonal dependence分离，产生多项相乘的式子（multiplicative seasonal ARIMA model）。其数学表达式如下：

1.4 Fourier Residual Modification

首先回忆一下傅里叶变换的概念：傅立叶变换，表示能将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。
此处傅里叶序列在FARIMA模型中是用于增加精度的，并且增加精度的作用很强;从本质上而言，傅里叶修正项的引入相当于也就是引入period的部分。
傅里叶修正项还可以用于处理多重周期：比如同时存在日周期和周周期的话（有几重周期就要有几个参数），就可以用AIC准则来获得所需要的周期参数。

以下博文介绍了如何利用fourier series来预测含有多重seasonality数据的方法和代码，可供参考：
https://content.pivotal.io/blog/forecasting-time-series-data-with-multiple-seasonal-periods

Fourier series确定的详细数学方法可参见参考资料中的pdf，这里给出一点概括；傅里叶修正项的数学表达式为：


结合其他网上的资料理解了以下，其中n-1代表T，是一个时间跨度中周期循环的次数（e.g. 一年平均有52.18个周）。F，即minimun deployment frequeny，代表着傅里叶修正项会有2F个三角函数式，F的计算由n得到；而三角函数的系数是由OLS估计得到的。
P.S. 虽然以上资料中给出的F是由n得到的，但是结合Rob J Hyndman的博文和其他资料后发现，傅里叶序列的阶数似乎并不一定等于(n-1)/2-1，而是遍历计算每一个阶数模型所对应的AIC值之后，选取AIC最小的阶数。或许可以理解为minimum deployment frequency是傅里叶函数阶数的一个特殊情况，但不一定是最有效的阶数取值方式。（个人理解，如有错误欢迎指正）

2 预测流程

之前一直对时间序列的预测方法流程的认知比较破碎，这次找到了一个较为完整的流程如下，个人觉得很有借鉴意义，后面的博文也大概按照这个流程图展开：

2.1 加载、检查数据

加载forecast和tseries程序包，注意两个都要加载，由于R自带的也有arima函数，因此如果没有加载的话会导致后面运行出的arima模型的结果值会有一些不同……

library(forecast)
library(tseries)

数据源自forecast package的作者Rob J Hyndman网站提供的汽油价格数据，这里导入时限制数据频率为周度是数据。
未来方便检测模型的预测结果，这里先删掉最后一年的数据，取得用作训练的时间序列数据tgas，然后用tsdisplay函观察数据特征。

gas <- ts(read.csv("https://robjhyndman.com/data/gasoline.csv",header=FALSE)[,1],freq=365.25/7,start=1991+31/365.25)
tgas<-ts(as.vector(gas[1:693]),frequency=365.25/7,start=1991+31/365.25)
plot(tgas)

从数据折线图的结果可以看出原来的gas数据存在明显的趋势和周期性波动（这样的周期性波动也反映在ACF图中），所以先作差去掉趋势，再检查数据是否平稳。

tgas_1<-diff(tgas)
#再次查看图形
tsdisplay(tgas_1)

从图上看，数据作差后基本围绕0在随机振动，基本平稳，进一步使用adf检验，看一下是否存在单位根。

adf.test(tgas_1)

#    Augmented Dickey-Fuller Test
#    data:  tgas_1
#    Dickey-Fuller = -11.924, Lag order = 8, p-value = 0.01
#    alternative hypothesis: stationary

#  adf结果显示拒绝单位根的假设


# 进一步观察一下数据大小分布，看看是不是正态分布的，结果基本满足，因而不再做进一步处理
hist(tgas_1,breaks=20)

此外，如果这一步中发现数据分布明显异常，也可以进行Box-Cox转换。
这里跳出来讲一下Box-Cox转换的定义：Box-Cox变换是统计建模中常用的一种数据变换，用于连续的相应变量不满足正态分布的情况。
当变量不满足正态分布时，对模型估计会有如下影响：
模型偏差和虚假交互：如果数据分布是不对称的（比如下图的分布，可见数据右端存在较多极值），这种不对称行为可能会导致回归模型的偏差。 因为数据较大时方差也较大，如此可能会导致虚假的相互作用，从而导致一个非常复杂的模型，其中包含许多（虚假和不真实的）变量间的交互作用。
Box-Cox转换的目的是让数据转为正态分布的，并且拥有稳定的方差。其数学表达式如下：

在这里Y是变换之前的数据，λ是一个待定变换参数，而Y(λ)则是变换之后的数据。对于不同的λ，所作的变换也不相同，所以Box-Cox变换是一族变换，它包括了平方根变换(λ=0.5)，对数变换(λ=0)和倒数变换(λ=-1)等常用变换。
为方便理解Box-Cox怎样把非正态分布的数据转为正态分布的数据，这里举λ=0，即对数变换的情况为例子：


（图片及实例来源：https://blog.minitab.com/blog/applying-statistics-in-quality-projects/how-could-you-benefit-from-a-box-cox-transformation）
当对数函数作用于数据时，较小数值之间的间隔会被扩大，而较大数据之间的间隔会被缩小，从而数据能转换为正态分布的形态。

2.2. 构建模型

因为好奇Seasonal ARIMA与FSARIMA的预测能力的差异，这里会同时计算两种模型的拟合结果。

2.2.1 Seasonal ARIMA

从原理上来讲，Seasonal ARIMA(p,d,q)(P,D,Q)[S]中季节性的阶数应该是可以通过观察ACF和PACF周期性变化的情况来确定的，但找了很多资料都觉得还是不太会看，所以这里还是直接使用auto.arima函数（该函数确定的结果可能是seasonal的，也可能是non-seasonal的）：

#ARIMA,参数通过auto.arima来确定，这样确定的结果可能是seasonal的，也可能是unseasonal的
autofit<-auto.arima(tgas)
autofit
#Series: tgas 
#ARIMA(0,1,2)(1,0,0)[52] 
#AIC=9909.91 

#储存auto.arima 生成的结果
fit2 <- arima(tgas,order=c(0,1,2),seasonal=list(order=c(1,0,0),period=52))

2.2.2 ARIMA with fourier modification

这里参考了Rob J Hyndman在博文中给出的代码，大致意思是构建一个循环，限定fourier modification的最高阶数，然后选出auto.arima()生成的函数与fourier modification结合之后AIC值最小的模型。
这里尝试了比较多次，发现auto.arima()生成的模型中，AIC最小的模型在加入了fourier modification得到的模型并不一定还是众模型中AIC最小的，因此先确定arima的阶数，再寻找最佳的fourier阶数不一定是最有效的方式。

bestfit_2 <- list(aicc=Inf)
for(i in 1:25)
{  
  fit <- auto.arima(tgas,xreg=fourier(tgas,K=i),seasonal=FALSE)
  if(fit$aicc <bestfit_2$aicc){
    bestfit<-fit  # arima(x = gas, order = c(0, 1, 2), xreg = fourier(gas, K = i))
    bestk<-i # 增加一行，以获得最佳的fourier order
    }
  else 
    break;
  # when there's more than one commands in one line, use ";"
}
bestfit # Regression with ARIMA(4,1,1) errors AIC=9710.95 
bestk #8
#储存循环生成的结果
fit3<-auto.arima(tgas,xreg=fourier(tgas,K=8),seasonal=FALSE)

2.3 模型效果评估

2.3.1 残差分析-自相关性

首先观察残差的自相关系数，当残差的自相关性较强时，这说明残差中还残留着可以用模型解释的部分，模型的效率不高，但并不是说模型的预测结果就是错误的。

acf(fit2$residual)

值得注意的是这里图像的横轴坐标有小数，这是因为acf（）返回的图像中横轴指代的是数据的一个周期，这里即52个数据，所以X=0.1就意味着5.2周；从图中可以看到acf的值基本在虚线范围内，说明seasonal ARIMA生成的模型中共线性问题并不严重。

acf(fit3$residual)

同上，残差自相关的问题通过检验。

2.3.2 残差分析-正态分布

观察残差是否是正态分布的。残差是否是正态分布对于预测来说并不是非常关键，但如果残差非正态分布，将导致预测的区间估计不那么准确。

checkresiduals(fit2)
#    Ljung-Box test

#data:  Residuals from ARIMA(0,1,2)(1,0,0)[52]
#Q* = 241.44, df = 101.36, p-value = 1.916e-13

#Model df: 3.   Total lags used: 104.357142857143

如上，checkresiduals()函数不仅返回residual的折线图、ACF图和直方图，还会返回ljung-box检验的结果；从折线图可以看出来残差没有明显的趋势，大致围绕0波动，但是根据acf图，残差的自相关系数呈现出周期波动的特点，这说明模型还有进一步优化的空间。此外，可以注意到checkresiduals()函数返回的acf图有横坐标为52.17，正好等于我们导入数据时定义的频率：365.25/7；因而该函数返回的acf图形中横坐标的一个单位代表一周；
再来看Ljung-Box test的检验结果，Ljung-Box test测试确定误差是否为IID（即白噪声），或是否背后有其他因素；误差或残差的自相关是否为非零。从本质上说，这是一个检验模型是否缺乏拟合能力的方法：如果残差的自相关非常小，我们就说模型没有显示出“显著的缺乏拟合力”。由于这里的H0是lack of fit,所以只要得到足够小的p值就可以证明模型的拟合能力。

上图为Ljung-Box Q值的表达式，其中T是时间序列的长度，rk代表残差的第k阶自相关系数，h是检验的滞后阶数。Q值服从自由度为h-k的卡方分布，Q值越大，残差序列自相关的特征就越强。
注意Ljung-Box检验中还需控制滞后阶数h，关于h的选择，Rob J Hydman建议在non-seasonal data中取阶数为10，取seasonal-data取阶数为seasonality period的两倍。如此操作是基于对检验效力的考虑。当h足够大时，才能够捕获任何有意义的或者引起麻烦的相关性。对于季节性数据，通常在残差中剩余的季节性滞后的倍数处具有相关性，因此我们希望至少包括两个季节周期的滞后项。
再倒回去看Seasonal ARIMA model Ljung-Box检验的结果，本次检验的阶数为104.36，恰为数据周期的两倍，而p值也远小于0.05，因此残差自相关性的检验基本通过。

接下来观察FARIMA模型checkresiduals的结果：

checkresiduals(fit3)

残差的折线图和直方图分布与Seasonal ARIMA大概相似，但acf图的表现更好，不像Seasonal ARIMA的ACF图显现出明显的周期性，acf的绝对值也更小。

2.4 预测

在网上找了一波资料，发现有用predict和forecast函数进行预测的两种方法，这里截取一下Rob J Hyndman对两种方法异同的解释：
二者给出的预测结果基本是一样的，但是forecast函数返回的是forecast函数， 相对于简单的list来说，画图，总结和分析会更方便一些。所以这边选用forecast函数。
先看Seasonal ARIMA的结果：

fc2 <- forecast(fit2,h=104)
plot(fc2)
lines(fc2$fitted,col="green")
lines(gas,col="red")

这里绿线是拟合值，红线是真实值，蓝线是预测值。预测时间段中浅蓝色和深蓝色的两个区域分别代表默认的80%以及95%置信区间内的数值边界。可以看到，利用历史数据“喂”出来的模型在第一年的预测结果中能捕捉到一定的数据周期变化的特征，但是到第二年的预测结果时就变得很平。
再来看FARIMA的结果：

fc3 <- forecast(fit3, xreg=fourier(tgas, K=8, h=104))
plot(fc3)
lines(fc3$fitted,col="green")
lines(gas,col="red")

在这张图中，可以看到FARIMA模型的预测结果能够表现更强的周期性，与最后一年实际数据的相差值也更小，预测更加精确。而且这样的周期性到第二年的预测结果中也仍旧保持着，并不像Seasonal ARIMA那样变得平缓。
这里再补充一下通过传统的绘制Acf, Pacf图，观察拖尾与截尾的情况得到的Non-Seasonal ARIMA模型的预测结果：

fit1<-arima(tgas,order=c(0,1,1))
fit1 <- auto.arima(tgas,seasonal=FALSE)
acf(fit1$residual)
acf(fit2$residual)
acf(fit2$residual[1:length(fit3$residual)])

Well ……这个预测结果真的是非常，非常地平。一开始真的很难接受在学校数学推导了大半个学期弄懂的模型就这玩意……怀疑自己弄错了，不过Non-seasonal ARIMA的假设就是无周期性、均值与方差有一定的稳定性，因此预测值也没有季节性的波动。由此可见Seasonality对于模型预测准确性的重要程度了。
季节性在商业中的应用很广，例如营业额、顾客流量、销售数目和unusual event study都是可能的例子，所以好好把这个数据特征和工具记心里叭。（第一次写博文，只是作为最近学习的一些记录。参考了很多资料以及代码，有些地方可能解释得不对，希望大家谅解&欢迎指正~如需转载，请注明出处）

参考资料：

最重要的代码参考：
RJH - Forecasting weekly data
https://robjhyndman.com/hyndsight/forecasting-weekly-data/
R语言 时间序列ARIMA模型方法
http://blog.csdn.net/gdyflxw/article/details/55509656

时间序列预测的流程、规范：
ARIMA modelling in R
（含时间序列预测的流程图，最后的实例也蛮不错，可以多看看）
https://robjhyndman.com/hyndsight/forecasting-weekly-data/

模型原理：
Seasonal ARIMA processes 数学原理
http://halweb.uc3m.es/esp/Personal/personas/amalonso/esp/TSAtema6.pdf
Seasonal ARIMA以及FARIMA的数学原理
（文章中写的是FSARIMA，因为本文设置的是Seasonal = False, 所以省去了S，自造了FARIMA的名字）
http://wseas.us/e-library/conferences/2013/Chania/AEBDb/AEBDb-08.pdf
Introduction to Forecasting with ARIMA in R
https://www.datascience.com/blog/introduction-to-forecasting-with-arima-in-r-learn-data-science-tutorials

一些细节解释：
利用fourier series来预测含有多重seasonality数据: Forecasting Time Series Data with Multiple Seasonal Periods
https://content.pivotal.io/blog/forecasting-time-series-data-with-multiple-seasonal-periods
box-cox变换 - 百度百科
https://baike.baidu.com/item/box-cox变换/10278422
box-cox变换的原理解释 - How Could You Benefit from a Box-Cox Transformation?
https://blog.minitab.com/blog/applying-statistics-in-quality-projects/how-could-you-benefit-from-a-box-cox-transformation
Evaluating the regression model
https://otexts.com/fpp2/regression-evaluation.html
Difference between forecast and predict function in R
https://stackoverflow.com/questions/31409591/difference-between-forecast-and-predict-function-in-r