Continuing from the previous time, we will use TOPIX historical data. In addition, I felt that there was data that seemed to have autocorrelation, so I prepared data on the number of foreign visitors to Japan every month. I used the one on the website of JTB Research Institute linked below. https://www.tourism.jp/tourism-database/stats/inbound/ At first glance, the data is seasonal, and the yen is depreciating in Abenomics, so the overall trend is upward.
wn = np.random.randn(50)
mu = 0.0
theta_1 = 1.0
map = mu + wn[1:] + theta_1 * wn[:-1]
The data was generated like this.
Considering the green line $ \ theta_1 = 0.0 $ as a reference, it is easy to understand the influence of the sign and value of $ \ theta_1
The main statistics are as follows.
mean :
Check the correlogram. In the plot above, the number of data is 50, but here it is created with 1,000 data.
Up to this point, it's an intuitive image. Only the first-order autocorrelation is significant, and when $ \ theta_1 <0 $, the autocorrelation is a negative value.
Next, let's look at a generalized version of the MA process.
The main statistics are as follows.
[Supplement] What is stationarity?
If less than $ \ quad $ holds, the process is weak stationarity.
Let's plot it. Here, $ y_0 = 0 $, $ c = 1 $, and $ \ sigma ^ 2 = 1 $.
AR1 = []
y = 0
c = 1
phi_1 = 0.5
for i in range(200):
AR1.append(y)
y = c + phi_1*y + np.random.randn()
The main statistics are as follows. However,
Next is the generalized AR process.
The main statistics are as follows.
The correlogram is as follows.
$ ARMA (p, q) $ has the following properties.
Let's actually plot $ ARMA (3,3) $. The parameters are as follows.
arma = [0,0,0]
wn = [0,0,0]
c = 1.5
phi_1 = -0.5
phi_2 = 0.5
phi_3 = 0.8
theta_1 = 0.8
theta_2 = 0.5
theta_3 = -0.5
for i in range(1000):
eps = np.random.randn()
y = c + phi_1*arma[-1] + phi_2*arma[-2] + phi_3*arma[-3] + eps + theta_1*wn[-1] + theta_2*wn[-2] + theta_3*wn[-3]
wn.append(eps)
arma.append(y)
Stationarity
(This is the same as the steady state condition of the AR process.)
AR characteristic equation :
Reversability
The MA process is reversible when the MA process can be rewritten into an AR (∞) process.
MA characteristic equation :
Least squares (OLS: ordinary least squares)
The OLS selects the parameters of the ARMA model to minimize the sum of squared residuals (SSR) between the actually observed values and the values estimated by the ARMA model.
You can solve the normal equations that are partially differentiated with respect to each parameter.
Maximum likelihood estimation (MLE)
The maximum likelihood method is the approach of selecting the parameter with the maximum log-likelihood.
From here, let's estimate the parameters for $ ARMA (3,3) $ actually created in 5. The statsmodels library is used because it is complicated to calculate by hand. See the URL below for details. https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html It seems that if you give the sample sample and the order of the model with order = (p, q) to the argument of ARMA, the parameter will be calculated by maximum likelihood estimation.
arma_model = sm.tsa.ARMA(arma, order=(3,3))
result = arma_model.fit()
print(result.summary())
ARMA Model Results
==============================================================================
Dep. Variable: y No. Observations: 1003
Model: ARMA(3, 3) Log Likelihood -1533.061
Method: css-mle S.D. of innovations 1.113
Date: Sun, 08 Dec 2019 AIC 3082.123
Time: 14:51:34 BIC 3121.409
Sample: 0 HQIC 3097.052
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 6.8732 0.410 16.773 0.000 6.070 7.676
ar.L1.y -0.4986 0.024 -21.149 0.000 -0.545 -0.452
ar.L2.y 0.5301 0.019 27.733 0.000 0.493 0.568
ar.L3.y 0.8256 0.020 41.115 0.000 0.786 0.865
ma.L1.y 0.7096 0.036 19.967 0.000 0.640 0.779
ma.L2.y 0.3955 0.039 10.176 0.000 0.319 0.472
ma.L3.y -0.4078 0.033 -12.267 0.000 -0.473 -0.343
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 1.0450 -0.0000j 1.0450 -0.0000
AR.2 -0.8436 -0.6689j 1.0766 -0.3933
AR.3 -0.8436 +0.6689j 1.0766 0.3933
MA.1 -0.6338 -0.8332j 1.0469 -0.3535
MA.2 -0.6338 +0.8332j 1.0469 0.3535
MA.3 2.2373 -0.0000j 2.2373 -0.0000
-----------------------------------------------------------------------------
And so on. We are able to estimate values that are reasonably close to the original parameters. A few notes, S.D. of innovations represents the dispersion of white noise, and to find the value of c from const,
print(result.params[0] * (1-result.arparams.sum()))
0.9824998883509097
Need to do as.
To select $ ARMA (p, q) $, first conservatively narrow down the $ p, q $ values from the autocorrelation and partial autocorrelation values, and then determine the optimal model based on the information criterion. take.
First, the k-th-order partial autocorrelation removes the influence of $ y_ {t-1}, \ cdots, y_ {t-k + 1} $ from $ y_t $ and $ y_ {tk} $. It is defined as the correlation of things.
Considering the above definition, the properties shown in the table below are clear.
model | Autocorrelation | 偏Autocorrelation |
---|---|---|
Decaying | ||
Decaying | ||
Decaying | Decaying |
This makes it possible to narrow down the model to some extent by looking at the structure of sample autocorrelation and sample partial autocorrelation.
Next is the information criterion (IC).
We have prepared data on TOPIX and the number of foreign visitors to Japan, but both are in the shape of not being steady. Therefore, I had no choice but to create data that excludes trends by subtracting the moving average for the 1996-2007 portion of the number of foreign visitors to Japan.
visit.head(3)
month | visits | |
---|---|---|
0 | 1996-01-01 | 276086 |
1 | 1996-02-01 | 283667 |
2 | 1996-03-01 | 310702 |
visit['trend'] = visit['visits'].rolling(window=13).mean().shift(-6)
visit['residual'] = visit['visits'] - visit['trend']
v = visit[visit['month']<'2008-01-01']
plt.plot(v['month'].values, v['visits'].values, label='original')
plt.plot(v['month'].values, v['trend'].values, label='trend')
plt.plot(v['month'].values, v['residual'].values, label='residual')
plt.legend()
plt.grid()
plt.title('Foreign Visitors')
plt.show()
As a trend, I decided to take a moving average of 6 months before and after. This is processed by rolling (window = 13) .mean () and shift (-6).
First, let's look at autocorrelation and partial autocorrelation. Very easy with the library.
fig = plt.figure(figsize=(8,5))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(v['residual'].dropna().values, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(v['residual'].dropna().values, lags=40, ax=ax2)
plt.tight_layout()
fig.show()
There is a strong seasonality that cannot be hidden. .. .. The correlation every 12 months is very strong. You can use statsmodels' arma_order_select_ic to estimate the optimal order of the ARMA model. For details, go to this link (https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.arma_order_select_ic.html)
sm.tsa.arma_order_select_ic(v['residual'].dropna().values, max_ar=4, max_ma=2, ic='aic')
Here, for $ ARMA (p, q) $, the AIC criterion is used with the maximum value of $ p $ being 4 and the maximum value of $ q $ being 2.
0 | 1 | 2 | |
---|---|---|---|
0 | 3368.221521 | 3363.291271 | 3350.327397 |
1 | 3365.779849 | 3348.257023 | 3331.293926 |
2 | 3365.724955 | 3349.083663 | 3328.831252 |
3 | 3361.660853 | 3347.156390 | 3329.447773 |
4 | 3332.100417 | 3290.984260 | 3292.800604 |
I feel that none of them change much, but the result is that $ p = 4, q = 1 $ is good.
arma_model = sm.tsa.ARMA(v['residual'].dropna().values, order=(4,1))
result = arma_model.fit()
plt.figure(figsize=(10,4))
plt.plot(v['residual'].dropna().values, label='residual')
plt.plot(result.fittedvalues, label='ARMA(4,1)')
plt.legend()
plt.grid()
plt.title('ARMA(4,1)')
plt.show()
It's like this. Personally, I felt that he was doing his best despite being a relatively simple model.
Recommended Posts