1 Introduction

Many companies issue stocks or shares that can be bought as an investment. A stock represents ownership on part of the company’s earnings and assets [6].

Like other investments, there is risk involved with buying stocks. The price of a stock fluctuates depending on how the company is doing. If the price increases after you bought it, you would profit, but if the price decreases, you would end up with a loss.

It would be reasonable to analyze past data in order to predict whether the stock price will increase or decrease in the short term. However, that is very difficult, if not impossible to do with complete accuracy. Instead, we study the volatility of the stock. The volatality of a stock measures the variance or standard deviation of the returns of the stock, which are the daily logged difference of prices in two consecutive days [7]. If the volatility is high, then the price is prone to fluctuating drastically. It is possible that the price would increase dramatically, but it is also possible that the price would plummet. As a result, a high volatility means that the stock has a high risk. On the other hand, if the volatility is low, the price would not fluctuate dramatically and would be more stable. A stock with low volatility is a less risky investment.

In this project, we will analyze the volatility of the Sony (SNE) stock. We will be fitting stochastic volatility models onto the demeaned returns of the stock and determine how well the model fits using the log-likelihood. We will create the model using the pomp package.

2 Description of the Dataset

The data used for this project is the daily adjusted closing price of the SNE stock from the past five years, starting from April 8th, 2013 to April 5th, 2018. The dataset is taken from Yahoo Finance [1].

We first read in the data. There are 1259 observations and 7 parameters. We will also take the log of the adjusted closing price.

Now, we plot the adjusted closing price of SNE versus the time and the log of the adjusted closing price versus time. Both plots show an upwards trend in the dataset. As time passes, the stock price of SNE steadily increases.

Next, we find the returns of the investment, where we take the difference between adjacent log prices. If \(y^*_{1:N}\) is the time series of weekly adjusted closing price, then the return time series \(z^*_{2:N}\) will be calculated using the equation

\[z^*_n=\log y^*_n - \log y^*_{n-1}\]

We plot the returns against time. The blue line on the plot represents the mean return value. We also find the demeaned returns of SNE, calculated by subtracting the mean from each value of \(z^*_n\).

3 Stochastic Volatility Model

We will first analyze the stochastic volatility model that we studied in class [4]. The model is written below.

\[ Y_n = exp(0.5H_n) \epsilon _n \]

\[ H_n = \mu_h (1-\phi) + \phi H_{n-1} + \beta _{n-1} R_n exp(-0.5H_{n-1}) + \omega _n \]

\[ G_n = G_{n-1} + \nu _n \]

where \(H_n\) is the log volatility, \(\beta _n = Y_n \sigma _n \sqrt{1-\phi ^2}\), \(\epsilon _n\) are iid \(N(0,1)\), \(\nu _n\) are iid \(N(0, \sigma _{\nu}^2)\), and \(\omega _n\) are iid \(N(0, \sigma _{\omega}^2)\).

3.1 Constructing the POMP model

We now build the POMP model of our stochastic volatility model. The code used is given from [4].

## Loading required package: pomp

We first declare the various names.

sne_statenames <- c("H","G","Y_state")
sne_rp_names <- c("sigma_nu","mu_h","phi","sigma_eta")
sne_ivp_names <- c("G_0","H_0")
sne_paramnames <- c(sne_rp_names,sne_ivp_names)
sne_covarnames <- "covaryt"

We now write the two rprocesses for our model. We also initialize the variables, with initial value of \(H\) being \(H_0\), initial value of \(G\) being \(G_0\) and initial value of Y_state being a random normal variable with mean 0 and variance \(\exp (0.5 H)\)

rproc1 <- "
  double beta,omega,nu;
  omega = rnorm(0,sigma_eta * sqrt( 1- phi*phi ) * sqrt(1-tanh(G)*tanh(G)));
  nu = rnorm(0, sigma_nu);
  G += nu;
  beta = Y_state * sigma_eta * sqrt( 1- phi*phi );
  H = mu_h*(1 - phi) + phi*H + beta * tanh( G ) * exp(-H/2) + omega;
"
rproc2.sim <- "
  Y_state = rnorm( 0,exp(H/2) );
 "
rproc2.filt <- "
  Y_state = covaryt;
 "
sne_rproc.sim <- paste(rproc1,rproc2.sim)
sne_rproc.filt <- paste(rproc1,rproc2.filt)

sne_initializer <- "
G = G_0;
H = H_0;
Y_state = rnorm( 0,exp(H/2) );
"
sne_rmeasure <- "
y=Y_state;
"
sne_dmeasure <- "
lik=dnorm(y,0,exp(H/2),give_log);
"

Now we would like to transform the parameters so that all are defined on the entire real line. We are given that \(0 < \phi < 1\), so we use a logistic transformation to \(\phi\) to define it in the real line. For both sigma parameters, we know that they must both be positive values, so we use log and exponential transformations for \(\sigma _\eta\) and \(\sigma _\nu\).

sne_toEstimationScale <- "
  Tsigma_eta = log(sigma_eta);
  Tsigma_nu = log(sigma_nu);
  Tphi = logit(phi);
"
sne_fromEstimationScale <- "
  Tsigma_eta = exp(sigma_eta);
  Tsigma_nu = exp(sigma_nu);
  Tphi = expit(phi);
"

We now build a pomp object that we will be using for filtering.

sne.filt <- pomp(data=data.frame(y=sne.ret.demeaned,
                                   time=1:length(sne.ret.demeaned)),
                   statenames=sne_statenames,
                   paramnames=sne_paramnames,
                   covarnames=sne_covarnames,
                   times="time",
                   t0=0,
                   covar=data.frame(covaryt=c(0,sne.ret.demeaned),
                                    time=0:length(sne.ret.demeaned)),
                   tcovar="time",
                   rmeasure=Csnippet(sne_rmeasure),
                   dmeasure=Csnippet(sne_dmeasure),
                   rprocess=discrete.time.sim(step.fun=Csnippet(sne_rproc.filt),delta.t=1),
                   initializer=Csnippet(sne_initializer),
                   toEstimationScale=Csnippet(sne_toEstimationScale), 
                   fromEstimationScale=Csnippet(sne_fromEstimationScale)
)

We will now simulate our stochastic volatility model using the simulate function. In addition, we build a pomp object for filtering for the simulations.

expit<-function(real){1/(1+exp(-real))}
logit<-function(p.arg){log(p.arg/(1-p.arg))}
params_test <- c(
  sigma_nu = exp(-4.5),  
  mu_h = -0.25,       
  phi = expit(4),     
  sigma_eta = exp(-0.07),
  G_0 = 0,
  H_0=0
)
sim1.sim <- pomp(sne.filt, 
                 statenames=sne_statenames,
                 paramnames=sne_paramnames,
                 covarnames=sne_covarnames,
                 rprocess=discrete.time.sim(step.fun=Csnippet(sne_rproc.sim),delta.t=1)
)
sim1.sim <- simulate(sim1.sim,seed=1,params=params_test)
sim1.filt <- pomp(sim1.sim, 
                  covar=data.frame(
                    covaryt=c(obs(sim1.sim),NA),
                    time=c(timezero(sim1.sim),time(sim1.sim))),
                  tcovar="time",
                  statenames=sne_statenames,
                  paramnames=sne_paramnames,
                  covarnames=sne_covarnames,
                  rprocess=discrete.time.sim(step.fun=Csnippet(sne_rproc.filt),delta.t=1)
)

After getting the simulated data, we plot it versus time. The plot is shown below. The plot in blue is the simulated demeaned returns, while the plot in black is the actual demeaned returns. We can see that the simulation might have overestimated the variance of the demeaned returns, so the model may not be the most accurate model to use for our dataset.

plot(Y_state~time, data=sim1.sim, type="l", col="blue", xlab="Time", ylab="Demeaned Returns", main="Simulated Demeaned Returns")
lines(zdem)

3.2 Filtering on the simulated data

In order to test the code, we specify different run levels for our program. For now, we run the program using run level 2.

run_level <- 2
sne_Np <-          c(100,1e3,2e3)
sne_Nmif <-        c(10, 100,200)
sne_Nreps_eval <-  c(4,  10,  20)
sne_Nreps_local <- c(10, 20, 20)
sne_Nreps_global <-c(10, 20, 100)

In order to run the program more efficiently, we do parallelization.

require(doParallel)
registerDoParallel()

Now, we will filter. We use the stew function to cache our results.

stew(file=sprintf("pf1-%d.rda",run_level),{
  t.pf1 <- system.time(
    pf1 <- foreach(i=1:sne_Nreps_eval[run_level], .export=c("sim1.filt","sne_Np","run_level"), 
                   .packages='pomp',
                   .options.multicore=list(set.seed=TRUE)) %dopar% try(
                     pfilter(sim1.filt,Np=sne_Np[run_level])
                   )
  )
},seed=493536993,kind="L'Ecuyer")
L.pf1 <- logmeanexp(sapply(pf1,logLik),se=TRUE)

3.3 Fitting the model to the SNE data

3.3.3 Local Search with Run Level 3

We now run the local search with run level 3, which should give more accurate results.

run_level <- 3
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3226    3232    3237    3243    3240    3278

The mean log-likelihood is 3243 and the maximum value is 3278. The values are slightly higher than those of the run level 2 local search, and are comparable with those of the run level 2 global search. The diagnostics plot and scatter plot also agree with the plots of the previous local search. Despite the better results than the local search of run level 2, the parameters do not converge as cleanly as the global search. For best results, doing a global search with run level 3 would be a potential extension.

4 Stochastic Volatility in Mean Model

We will now build another stochastic model for our demeaned returns dataset. This time, we will use the stochastic volatility in mean model.

According to [5], the model can be defined as follows.

\[ Y_n = d \sigma ^2 \exp (H_n) + \sigma \exp (0.5 H_n) \epsilon _n \]

\[ H_n = \phi H_{n-1} + \sigma _\nu \nu _n\]

Where \(Y_n\) are the demeaned returns, \(\epsilon _n\) and \(\nu _n\) are iid \(N(0,1)\) random variables, \(H_n\) is the volatility, and \(0 < \phi < 1\).

4.1 Constructing the POMP model

Like for the first model, we will write down the variable names

sne_statenames <- c("H","Y_state")
sne_rp_names <- c("d", "phi", "sigma_eta", "sigma")
sne_ivp_names <- c("H_0")
sne_paramnames <- c(sne_rp_names,sne_ivp_names)
sne_covarnames <- "covaryt"

We write down the two rprocesses for the model. For the Y_state, we write it as a random normal variable of mean \(d \sigma ^2 \exp (H_n)\) and variance \(\sigma \exp (0.5 H_n)\) [8]. This is because \(\sigma \exp (0.5 H_n)\) is multiplied by \(\epsilon_n\), which is a random normal variable of mean 0 and variance 1. The second term of \(Y_n\) affects the variance, while the first term would be the mean.

rproc1 <- "
double eta;
eta = rnorm(0, 1);
H = phi * H + sigma_eta * eta;
"
rproc2.sim <- "
Y_state = rnorm(d * sigma * sigma * exp(H), sigma * exp(0.5*H));
"
rproc2.filt <- "
Y_state = covaryt;
"
sne_rproc.sim <- paste(rproc1,rproc2.sim)
sne_rproc.filt <- paste(rproc1,rproc2.filt)

sne_initializer <- "
H = H_0;
Y_state = rnorm(d * sigma * sigma * exp(H), sigma * exp(0.5*H));
"
sne_rmeasure <- "
y=Y_state;
"
sne_dmeasure <- "
lik=dnorm(y, d * sigma * sigma * exp(H), sigma * exp(0.5*H),give_log);
"

We also transform the parameters, like in the previous model building. Since \(\sigma _\eta\) and \(\sigma\) are postive parameters, we use a log and exponent transform. Since \(0 < \phi < 1\), we use a logistic transformation.

sne_toEstimationScale <- "
Tsigma_eta = log(sigma_eta);
Tsigma = log(sigma);
Tphi = logit(phi);
"
sne_fromEstimationScale <- "
Tsigma_eta = exp(sigma_eta);
Tsigma = exp(sigma);
Tphi = expit(phi);
"

We now give initial values for the parameters.

params_test <- c(
  d = .5,  
  phi = 0.99,  
  sigma = 0.99,
  sigma_eta = 0.03,
  H_0=0
)

After the initializing of parameters, the code will be the same as the code of the previous model. We do a local search and give the summary of the log-likelihoods found.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3256    3277    3278    3277    3282    3285

We see that the mean log-likelihood is 3277 and the max value found is 3285. Now we give the scatter plot and diagnostics plot.

Here, we see that the value of \(\sigma\) is very close to 0. This is problematic since the demeaned returns is defined by \(Y_n = d \sigma ^2 \exp (H_n) + \sigma \exp (0.5 H_n) \epsilon _n\). If \(\sigma\) is close to 0, we would get \(Y_n=0\) as well, resulting in a value of 0 for the demaned returns. The log-likelihood is high in this case because the mean of the demeaned returns is 0.

4.2 Modified POMP Model

To help solve the problem, we would fix the value of \(\sigma\) to be 1. Then the stochastic volatility in mean model is defined by

\[ Y_n = d \exp (H_n) + \exp (0.5 H_n) \epsilon _n \]

\[ H_n = \phi H_{n-1} + \sigma _\nu \nu _n\]

We now declare the variable names and rprocesses.

sne_statenames <- c("H","Y_state")
sne_rp_names <- c("d", "phi", "sigma_eta")
sne_ivp_names <- c("H_0")
sne_paramnames <- c(sne_rp_names,sne_ivp_names)
sne_covarnames <- "covaryt"

rproc1 <- "
double eta;
eta = rnorm(0, 1);
H = phi * H + sigma_eta * eta;
"
rproc2.sim <- "
Y_state = rnorm(d * exp(H), exp(0.5*H));
"
rproc2.filt <- "
Y_state = covaryt;
"
sne_rproc.sim <- paste(rproc1,rproc2.sim)
sne_rproc.filt <- paste(rproc1,rproc2.filt)

sne_initializer <- "
H = H_0;
Y_state = rnorm(d * exp(H), exp(0.5*H));
"
sne_rmeasure <- "
y=Y_state;
"
sne_dmeasure <- "
lik=dnorm(y, d * exp(H), exp(0.5*H),give_log);
"
sne_toEstimationScale <- "
Tsigma_eta = log(sigma_eta);
Tphi = logit(phi);
"
sne_fromEstimationScale <- "
Tsigma_eta = exp(sigma_eta);
Tphi = expit(phi);
"

We filter and simulate. We also plot the simulation, comparing it to the actual dataset.

params_test <- c(
  d = .5,  
  phi = 0.99,     
  sigma_eta = 0.03,
  H_0=0
)

We can see that the simulation overestimates the values of the demeaned returns.

4.2.1 Local Search

We continue on by performing the local search.

##    user  system elapsed 
##    0.15    0.04  982.64
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3223    3226    3229    3228    3231    3234

The program ran for 982.64 seconds, or 16.377 minutes. We obtain a mean log-likelihood value of 3228 and a max value of 3234. We now plot the scatterplot and diagnostics.

Looking at the scatterplot, we see that \(\phi\) tends to go towards a value near 1. From the diagnostics, we see that there are times where the effective sample size drops significantly, when the dataset at that time oscillates more dramatically. We see the log-likelihood quickly converges to a value of about 3200, \(\phi\) converges near 1, and \(\sigma _\eta\) converges around 0.3.

4.2.2 Global Search

We now perform the global search. We give a range of values for the parameters, and run the program below.

##    user  system elapsed 
##    0.19    0.21  964.27
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3216    3222    3228    3228    3232    3238

The program ran for 964.27 seconds, or 16.07 minutes. The mean log-likelihood found is 3228 and the maximum value found is 3238.

Looking at the scatterplot, we see that \(\phi\) tends to go near a value of 1, and \(\sigma _\eta\) being around 0.3. We also see a linear relationship between the log-likelihood and \(H_0\). As \(H_0\) increases, the log-likelihood decreases.

Looking at the diagnostics plot, we see similar results for the effective sample size as for the local search. The log-likelihood converges quickly to a little higher than 3200, while \(\phi\) converges near 1 and \(\sigma _\eta\) converges around 0.3.

4.2.3 Local Search with Run Level 3

Like for the previous model, we run the local search with run level 3.

run_level <- 3
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3224    3230    3232    3232    3234    3238

The local search results in a mean log-likelihood value of 3232 and a maximum value of 3238. According to the scatterplot and diagnostic plot, we see that \(\phi\) converges to a value of 1, \(\sigma _\eta\) converges near 0.3, and the log-likelihood values are consistent. The results are improved from the previous local search and slightly improved in terms of mean value compared to the global search.

5 Benchmark Likelihood from GARCH model

Having obtained log-likelihood values for our models, we will compare them to that of a simpler model that can be used for our dataset. In this case, we will use the GARCH(1,1) model.

The GARCH(1,1) model has the form

\[ Y_n = \epsilon _n \sqrt{V_n} \] \[ V_n = \alpha _0 + \alpha _1 Y_{n-1}^2 + \beta _1 V_{n-1} \]

We will find the log likelihood of the GARCH(1,1) model.

require(tseries)
fit.garch.benchmark <- garch(zdem,grad = "numerical", trace = FALSE)
L.garch.benchmark <- logLik(fit.garch.benchmark)
L.garch.benchmark
## 'log Lik.' 3160.049 (df=3)

The log-likelihood obtained from the GARCH(1,1) model is 3160.049 with 3 degrees of freedom.

The max log-likelihood obtained from the first stochastic volatility model is 3278 while the max log-likelihood obtained from the stochastic volatility in mean model is 3238. Based on the value of the log-likelihood, we can conclude that our models perform better than the GARCH(1,1) model.

6 Conclusion

In this project, we studied the volatility demeaned returns of the Sony stock from the past five years. In order to do this, we fitted two kinds of volatility models to our dataset and compared the log-likelihood values. We first created the pomp model and ran local and global searches for each model.

As a result, we obtained a maximum log-likelihood value of 3278 from the first stocastic volatility model, and a value of 3238 from the stochastic volatility in mean model. After obtaining the log-likelihood values, we compared them to that of the GARCH(1,1) model in order to see how well how models performed compared to a simpler but easily coded model. The GARCH(1,1) resulted a log-likelihood of 3160.049, which is lower than the log-likelihood found with our models. We can conclude that our models performed better than the bench value.

6.1 Extensions

In this project, we set the run level to 2 for most of the searches, while we only used run level 3 for the local searchs of the two volatility models. A further analysis can start with using the run level 3 on global searches and compare the results with the results described in this project. In addition, with further computing power, we may also have a higher run level with larger number of particles and iterations.

7 Sources

[1] https://finance.yahoo.com/quote/SNE/history?period1=1364702400&period2=1522468800&interval=1wk&filter=history&frequency=1wk

[2] https://ionides.github.io/531w18/

[3] http://ionides.github.io/531w16/final_project/index.html

[4] https://ionides.github.io/531w18/14/notes14.html

[5] Koopman, Siem Jan and Uspensky, Eugenie Hol. “The Stochastic Volatility in Mean model:Empirical evidence from international stock markets.”

[6] https://www.investopedia.com/terms/s/stock.asp

[7] https://www.investopedia.com/terms/v/volatility.asp

[8] https://ionides.github.io/531w16/final_project/Project08/Final/Stats_531_Final_Project.html