R语言 - seasonal ARIMA与带傅里叶修正项的ARIMA预测及比较

最近接触到了带傅里叶修正项的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