1. Introduction

The scientific question motivating my work is how to predict the stock price using time series, and based on our prediction models, how the stock price will behave in future. However, the common view about stock price prediction is that the daily returns are independent so that there is no chance for someone to make profit with a precise prediction model. The underlying reason is that if the prediction model tells that the stock price will go up and then people will buy more which will push up the price and eliminate the favorable situation (Lecture notes 15, Stats 531, Ed Ionides). Instead, in this project, we will study the volatility of the returns which might be predictable. Based on volatility of major companies’ stock price, we aim to fit applicable time series models which may predict the future behavior of the stock price.

2. Data

In this project, we use the adjusted close price for JP Morgan Chase & Co. from 2000 to 2016 from www.finance.Yahoo.com. The goal is to identify the best model to make a prediction for the volatility of this price. The reason we use he adjusted close price is that it is more meaningful when examining historical returns “because it gives analysts an accurate representation of the firm’s equity value beyond the simple market price and it accounts for all corporate actions such as stock splits, dividends/distributions and rights offerings” (www.investopedia.com). We will use the data from 2000 to 2013 as the training data to build the model and use the rest as test data to see that whether the model gives a reasonable prediction.

Firstly, we have a glimpse of the raw data of adjusted close price of JP Morgan. From the plot of training data, we can see that overall there is an increasing trend with some obvious fluctuation around the year of 2002 and 2008 which due to the financial crisis in Argentina and America, respectively. Additionally, the histogram shows that the distribution of adjusted close price is approximately normal.

Next, we calculate the difference of the log adjusted close price. As shown below, the large log differences occur around the year of 2002 and 2008 and the mean of the log differences is about 0. Moreover, according to the QQ-plot, our data roughly follows a normal distribution with slight long right tail.

3.Fit Models

3.1 ARIMA Models

We first take the advantage of ARIMA model, but should expect a terrible fit due to the common view about the independence of returns while ARIMA model assume there is correlation between daily returns although here we are dealing with the volatility. We do the model selection based on BIC values. According to the BIC table below, we choose ARIMA(1,1,1).

MA0 MA1 MA2 MA3 MA4 MA5
AR0 -15276.86 -15295.55 -15291.09 -15285.63 -15278.34 -15273.84
AR1 -15293.76 -15295.75 -15288.82 -15280.44 -15272.65 -15266.42
AR2 -15289.81 -15288.82 -15280.65 -15272.59 -15264.34 -15284.38
AR3 -15284.36 -15280.73 -15272.66 -15264.45 -15256.43 -15277.03
AR4 -15276.95 -15272.54 -15264.10 -15256.65 -15253.12 -15267.52
AR5 -15272.94 -15266.65 -15284.95 -15276.90 -15268.67 -15261.02
## [1] -15295.75

After fitting the ARIMA\((1,1,1)\) model, we do the acf and pacf plots, the correlation is not significant as shown for most lags, so we can confirm that the model has no seasonal period. There is also a significant correlation at lag1 in acf. The pacf shrinks slowly. These information suggests that a simple model without seasonality might be adequate. Thus, we can specify our multiplicative ARIMA(1,1,1). The log likelihood for model m1 is 7652.29.

3.2 ARIMA(1,1,1) Diagnostic check

The plot of residuals suggests some major irregularities with the model (relatively large residuals for the years of 2002 and 2008 again). The histogram of residuals looks normal. Acf and pacf show that residuals of the fitted model have some serial correlation, and thus we cannot accept that they are independent. This undesirable fitting might be due to the violation of underlying assumption for ARIMA models. To get a better fitting, we may consider to use other models.

3.3 GARCH Models

Our next approach is using GARCH models. We try to fit the residuals of fitted ARIMA model into a GARCH model and see the performance of residuals of GARCH. The eacf of absolute residuals and squared residuals below suggest GARCH(1,1) and GARCH(1,2) respectively. The AIC for GARCH(1,1) is -17699.21 and AIC for GARCH(1,2) is -17657.67. Here we choose GARCH(1,1) which has the smaller AIC and also GARCH(1,1) model is a popular choice (Cowpertwait and Metcalfe 2009, Edward Ionides).

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 x o o o o o o o o o o  o  o  x 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x x x o o o o o o o o  o  o  o 
## 4 x x x x o o o o o o o  o  o  o 
## 5 x x x x x o o o o o o  o  o  o 
## 6 x x x o x x o o o o o  o  o  o 
## 7 x x x o o x x o o o o  o  o  o
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 x x o o x x o x x x o  x  o  x 
## 2 x x x o o x x x x x x  x  o  o 
## 3 x x o o o x x o o x o  o  x  o 
## 4 x x x o x x o o o x o  o  x  o 
## 5 x x x x x x o o o x o  o  x  o 
## 6 x x x x o x x o o o o  o  o  o 
## 7 x x x x o o x x o o o  o  o  o
## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     6.782312e-04     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -1.107e+04
##      1    6 -1.109e+04  2.14e-03  3.38e-03  1.0e-03  3.8e+09  1.0e-04  6.35e+06
##      2    7 -1.109e+04  6.05e-05  8.94e-05  9.7e-04  2.0e+00  1.0e-04  2.88e+02
##      3   14 -1.124e+04  1.32e-02  2.29e-02  5.0e-01  2.0e+00  1.0e-01  2.79e+02
##      4   21 -1.125e+04  1.02e-03  4.18e-03  3.4e-04  4.6e+00  1.0e-04  2.08e+01
##      5   22 -1.126e+04  3.50e-04  2.70e-04  3.2e-04  2.0e+00  1.0e-04  2.56e+01
##      6   23 -1.126e+04  3.08e-05  3.51e-05  3.3e-04  2.0e+00  1.0e-04  2.77e+01
##      7   30 -1.134e+04  7.43e-03  7.02e-03  3.2e-01  2.0e+00  1.1e-01  2.72e+01
##      8   32 -1.152e+04  1.52e-02  1.57e-02  3.9e-01  2.0e+00  2.1e-01  2.17e+03
##      9   34 -1.180e+04  2.35e-02  1.95e-02  2.2e-01  1.9e+00  2.1e-01  1.22e+00
##     10   36 -1.185e+04  4.52e-03  5.41e-03  3.5e-02  2.0e+00  4.2e-02  5.30e+01
##     11   38 -1.191e+04  5.03e-03  5.07e-03  3.2e-02  2.0e+00  4.2e-02  6.08e+01
##     12   40 -1.196e+04  4.09e-03  4.39e-03  3.0e-02  2.0e+00  4.2e-02  3.41e+01
##     13   47 -1.196e+04  3.93e-06  8.59e-06  4.0e-07  2.7e+01  5.6e-07  4.08e-01
##     14   48 -1.196e+04  8.53e-08  8.60e-08  3.9e-07  2.0e+00  5.6e-07  5.25e-01
##     15   57 -1.199e+04  2.83e-03  6.21e-03  4.9e-02  2.0e+00  7.4e-02  5.24e-01
##     16   58 -1.202e+04  2.58e-03  3.10e-03  3.8e-02  9.8e-01  7.4e-02  3.30e-03
##     17   60 -1.205e+04  1.90e-03  1.72e-03  3.7e-02  0.0e+00  7.4e-02  3.38e-03
##     18   62 -1.205e+04  6.06e-04  7.37e-04  6.6e-03  1.9e+00  1.5e-02  1.25e-02
##     19   65 -1.206e+04  6.57e-04  7.82e-04  1.7e-02  1.4e+00  3.9e-02  6.83e-03
##     20   66 -1.207e+04  5.84e-04  7.72e-04  1.7e-02  8.7e-01  3.9e-02  1.41e-03
##     21   72 -1.207e+04  8.44e-06  2.93e-05  1.3e-07  8.3e+00  2.4e-07  1.31e-04
##     22   73 -1.207e+04  2.15e-06  2.05e-06  9.6e-08  2.0e+00  2.4e-07  3.89e-05
##     23   74 -1.207e+04  5.01e-09  5.08e-09  1.0e-07  2.0e+00  2.4e-07  3.80e-05
##     24   83 -1.207e+04  1.75e-05  2.80e-05  3.3e-03  6.4e-01  7.8e-03  3.79e-05
##     25   84 -1.207e+04  1.36e-07  6.73e-07  1.6e-04  0.0e+00  3.2e-04  6.73e-07
##     26   85 -1.207e+04  2.42e-07  2.51e-07  1.5e-04  4.3e-01  3.2e-04  2.84e-07
##     27   86 -1.207e+04  2.72e-08  2.86e-08  2.3e-05  0.0e+00  4.9e-05  2.86e-08
##     28   87 -1.207e+04  3.26e-10  4.11e-10  1.4e-05  0.0e+00  3.4e-05  4.11e-10
##     29   88 -1.207e+04  1.54e-11  1.33e-11  2.3e-06  0.0e+00  5.9e-06  1.33e-11
## 
##  ***** RELATIVE FUNCTION CONVERGENCE *****
## 
##  FUNCTION    -1.206889e+04   RELDX        2.317e-06
##  FUNC. EVALS      88         GRAD. EVALS      30
##  PRELDF       1.325e-11      NPRELDF      1.325e-11
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    1.483531e-06     1.000e+00    -4.646e+00
##      2    7.257394e-02     1.000e+00     1.155e-03
##      3    9.279388e-01     1.000e+00     1.125e-03

3.4 GARCH(1,1) Diagnostic check

The histogram of residuals of GARCH(1,1) shows perfect normality. The qq plot of residuals of GARCH(1,1) still suggests some heavily-tailed distribution, but it slightly improved compared with the qq plot of reiduals of fitted ARIMA model. The acf plots of absolute and squared residuals of GARCH(1,1) show that both absolute residuals and squared residuals have little serial correlation. We may conclude that GARCH(1,1) gives a better fit than ARIMA(1,1,1)

All p-values of generalized portmanteau tests are higher than 5%, which confirms that the standardized residuals from the fitted GARCH(1,1) model are uncorrelated. Hence, they may be independent. Also, We get the periodogram of the residuals of the GARCH(1,1) model shown below. The periodogram of the residuals of the GARCH(1,1) has no significant trend, which suggests white noise. This 3-parameter model has a maximized log likelihood of -36373.6.

3.5 POMP Models

The last method we use is POMP model. Since the plot of volatility is quite similar to the plot in the lecture notes 15(Edward Ionides), we first try the model as following:

\[ Y_n = \exp\{H_{n}/2\} \epsilon_n,\\ H_n = \mu_h(1-\phi) + \phi H_{n-1} + \beta_{n-1}R_n\exp\{-H_{n-1}/2\} + \omega_n,\\ G_n = G_{n-1}+\nu_n,\\ \] where \(\beta_n=Y_n\sigma_\eta\sqrt{1-\phi^2}\), \(\{\epsilon_n\}\) is an iid N(0,1) sequence, \(\{\nu_n\}\) is an iid N(0,\(\sigma_{\nu}^2\)) sequence and \(\{\omega_n\}\) is an iid N(0,\(\sigma_\omega^2\)) sequence.

Here we choose the algorithmic parameters based on the MIF2 convergence diagnostics. However, after trying several seemly reasonable parameters, we fail to reach a convincing result. Under this circumstance, we decide to fix some of those algorithmic parameters and see whether other parameters and log likelihood converge. Here, we let \(G_0=0,H_0=0,\mu_h=-8,\sigma_\eta=1\). We deduce those information based on a level 3 diagnostics and the reason we fix these four parameters is that they converge fast before the other two parameters converge. Since \(\phi\) and \(\sigma_\nu\) carry the most important information in our financial model, it is acceptable to only tune these two parameters. Moreover, we change the cooling.fraction.50 to 0.1 to make parameters converge faster. The MIF2 convergence diagnostics plots below indicate that all log likelihood, \(\phi\) and \(\sigma_\nu\) converges after about 10 MIF iteration.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  8875.2  8883.3  8885.7  8884.8  8888.0  8890.1

The profile likelihood of \(\sigma_\nu\) is shown below. From the plot, we find that the log likelihood achieves its maximum when \(\sigma_\nu\) is between 0.001 and 0.002, which is really close to zero. This might be due to those fixed algorithmic parameters and the searching area (from 0 to 0.2) but the result is still acceptable since the 95% confidence interval for \(\sigma_\nu\) does not contain zero.

We apply the same method to \(\phi\) and the profile likelihood is shown below. Again, the maximal log likelihood is achieved around the boundary of \(\phi\) when the searching area is from 0.8 to 0.99. Although the result is not perfectly desirable, we conclude that the 95% confidence interval for \(\phi\) is [0.97, 0.98].

4. Simulation using POMP Model

After obtaining the estimates of the algorithmic parameters, we can evaluate our model by comparing the original data with simulation results given by our model. If all simulation plots follow the patterns of the original data, we can say that our model give a reasonable fit and will provide a meaningful predication. As shown below, the blue plot stands for our original data and three black plots are the results of our simulation. Generally, the patterns of the simulation results are similar to the pattern of our original data. Moreover, we plot the confidence intervals (red in the plots) for the original data (blue in the plots) with simulated volatility. As we can see, most data are contained in our confidence interval except the data from the year of 2003, which are identified as a year with financial crisis. In all, the model provides a acceptable fit with some imperfections.

5.Conclusion and future work

Comparing all three models above, we may conclude that POMP model is the best one. ARIMA(1,1,1) is undesirable is due to the fact that it assume that there is a correlation between consecutive data. GARCH(1,1) seems to be better than ARIMA(1,1,1) but it is hard to be interpreted and there is no useful information we can draw from it. By constructing a POMP model, we can get the estimated parameters, such as \(Y_n\) and \(\sigma_\nu\), which are meaningful in our financial model and provide detailed explanation for the volatility. The advantage of POMP model is that we can still make improvement of the model by tuning the parameters or enlarge the number of iteration. Due to the considerable amount of computation and limited time, we are not able to identify the best estimates of all algorithm parameters (which are the things we can do next). In a word, we will stay with the POMP model and try more combinations of algorithm parameters to find the best ones to interpret our data.

6.Reference

[1]Lecture notes/Homework Solution of Stats 531 (Winter 2016) ‘Analysis of Time Series’, instructor: Edward L. Ionides (http://ionides.github.io/531w16/)

[2]Investopedia (http://www.investopedia.com/)

[3]Bretó, C. 2014. On idiosyncratic stochasticity of financial leverage effects. Statistics & Probability Letters 91:20–26.