# Volatility models applied to geophysics and high frequency financial market
data

Maria C. Mariani111Department of Mathematical Sciences and Computational
Science Program, University of Texas at El Paso. Md Al Masum
Bhuiyan222Computational Science Program, University of Texas at El Paso. Osei
K. Tweneboah333Computational Science Program, University of Texas at El Paso
Hector Gonzalez-Huizar444Department of Geological Sciences, University of
Texas at El Paso. and Ionut Florescu555School of Business, Stevens Institute
of Technology, Hoboken.

###### Abstract

This work is devoted to the study of modeling geophysical and financial time
series. A class of volatility models with time-varying parameters is presented
to forecast the volatility of time series in a stationary environment. The
modeling of stationary time series with consistent properties facilitates
prediction with much certainty. Using the GARCH model for financial stationary
data and the stochastic volatility model for geophysical stationary data, we
forecast one-step-ahead suggested volatility with $\pm 2$ standard prediction
errors, which is enacted via Maximum Likelihood Estimation. We compare the
stochastic volatility model relying on the filtering technique as used in the
conditional volatility with the GARCH model. We conclude that the stochastic
volatility is a better forecasting tool than GARCH (1,1), since it is less
conditioned by autoregressive past information.

Keywords: ADF test; KPSS test; Financial time series; Geophysical time series;
GARCH model; Maximum likelihood estimation; Stochastic volatility model;
Seismogram

##  1 Introduction

Forecasting of time series with estimation of time-varying parameters is very
important in modeling the dynamic evolution of the volatility. It is assumed
that a model that attracts the attention of investors can potentially be used
to predict key variables, for instance, returns, volatility, and volume of
stock market. It is to be noted that the development of forecasting
methodologies in geophysics helps us to identify the type of source that
generates a recorded seismic signal. This type of methodologies is generally
applied to various fields, such as finance, geophysics, and safety of power
system [1]. So, a reliable technique of forecasting, including the related
time information, is essential to construct less risky portfolios or to make
higher profits.

Financial time series manifests typical non-linear characteristics, and they
involve volatility clustering where the returns indicate their dynamism. In
this study, we develop a volatility forecasting method in which the logarithm
of the conditional volatility follows an autoregressive time series model.
R.F. Engle’s paper [2] introduced the autoregressive conditional
heteroskedasticity (ARCH) model to express the conditional variance of
available returns as a function of previous observations. Few years later, S.
Bollerslev [3] modified this concept and generalized the ARCH (GARCH) model
that allows the conditional variance to depend on the previous conditional
variance as well as on the squares of previous returns. In other words, the
system volatility in GARCH model is driven by the observed values in a pre-
deterministic fashion. In fact, over the past few decades, a considerable
amount of deterministic models has been suggested to forecast the volatility
in finance. The reason is that they are very simple, and help to account for
clustered errors and nonlinearity issues. In the present study, we first
propose a continuous-time stationary and GARCH (1,1) process that is useful in
the analysis of high frequency financial time series.

It is now widely believed that the measurements of a sequence of geophysics
and finance are stochastically dependent on the time needed. In other words,
there is a correlation among the numbers of data points at successive time
intervals. In Ref. [5] and [6] , the authors used stochastic models to
describe a unique type of measurement dependence in geophysical and financial
data. It has been observed that the data may follow different behaviors over
time, for instance, the mean reversion and fluctuation of power spectrum. Such
observations are unlike those of the classical modeling foundations. But the
concept of time-dependent observations suggests that the current information
needs to be evaluated on the basis of its past behavior [4]. This behavior of
time series makes it possible to effectively forecast volatility and to obtain
some stylized facts, namely, time-varying volatility, persistence, and
clustering.

A distinctive feature of the time series is that the deterministic model i.e.
GARCH does not allow for a full statistical description of volatility [7].
When there are high fluctuations in the time series, the GARCH (1,1) model
predicts the volatility arbitrarily since it cannot capture the high volatile
nature of the data. So in this work, we propose a stochastic model and
filtering technique as a way to estimate parameters. We therefore study a
continuous-time stationary sequences corresponding to seismograms of mining
explosions and small intraplate earthquakes to forecast the volatility by
using estimated parameters. These stationary sequences are very effective to
capture the characteristic of time-varying parameters in an appropriate way
[8]. We determined the adequacy and stationarity of the data by computing the
estimated standard error and some powerful tests respectively, which will be
discussed in Sections 4 and 5 of this paper. The main difficulty of SV model
is to fit it into the data (with higher accuracy in a stochastic process),
since their likelihood estimations involve numerical integration over higher
dimensional intractable integrals [9], whose maximization is rather complex.

The paper is organized as follows: Section 2 describes the methodology of
forecasting the volatility of time series. The techniques for estimating the
model parameters will also discussed. Section 3 is devoted to the background
of data. We discuss some essential parts of the financial and geophysical data
properties. In section 4 we perform tests that analyzes the stationarity of
the time series. Sections 5 and 6 provide the results when our models are
applied to the data sets. These sections also study the suitability of our
model with reference to the estimation of model parameters and the prediction
of the volatility of time series. Finally, section 7 contains the conclusion
and makes a comparison between the models.

##  2 Methodology

This section describes the volatility models that will be used to forecast the
time series regarding finance and geophysics. We will discuss some techniques
that will be convenient to estimate the parameters of the proposed models.

###  2.1 Filtering Approach

The state space model is defined by a relation between the $m$-dimensional
observed time series, ${\bf y}_{t}$, and the $n$-dimensional state vector
(possibly unobserved), ${\bf x}_{t}$ [10]. An observed equation is driven by
the stochastic process as follows:

${\bf y}_{t}={\bf Z}_{t}{\bf x}_{t}+{\bf\epsilon}_{t},$ (1)

where ${\bf Z}_{t}$ is a $m\times n$ observation matrix, ${\bf x}_{t}$ is a
vector of $n\times 1$, and ${\bf\epsilon}_{t}$ is a Gaussian error term
(${\bf\epsilon}_{t}\sim N(0,{\bf\beta}_{t})).$
The unobservable vector ${\bf x}_{t}$ is generated from the transition
equation which is defined as:

${\bf x}_{t}={\bf Tx}_{t-1}+{\bf\delta}_{t},$ (2)

where ${\bf T}$ is a $n\times n$ transition matrix and ${\bf\delta}_{t}\sim$
i.i.d $N(0,{\bf\zeta_{t})}$. We assume that the process starts with a Normal
vector $\bf{x_{0}}$. From Eqs. (1) and (2), we make estimation for the
underlying unobserved data ${\bf x}_{t}$ from the given data ${\bf
Y}_{m}=\\{y_{1},\ldots,y_{m}\\}$. When $m=t$, the process is called filtering.

###  2.2 Likelihood Approximation

Let $\varphi$ denote the parameters of the state space model, which are
embedded in the system matrices ${\bf Z}_{t},\bf{T},{\bf\beta}_{t},$ and
${\bf\zeta}_{t}$. These parameters are typically unknown, but estimated from
the data $Y={y_{1},\ldots,y_{m}}$.

The likelihood $L(\varphi\big{|}Y)$ is a function that assigns a value to each
point in the parameter space $\Delta$ which suggests the likelihood of each
value in generating the data. However, the likelihood is proportional to the
joint probability distribution of the data as a function of the unknown
parameters. The maximum likelihood estimation means the estimation of the
value of $\varphi\in\Delta$ that is most likely to generate the vector of the
observed data $y_{t}$ [11]. We may represent this as:

$\hat{\varphi}_{MLE}=\smash{\displaystyle\max_{\varphi\in\Delta}}L(\varphi\big{|}Y)=\smash{\displaystyle\max_{\varphi\in\Delta}}L_{Y}(\varphi)=\smash{\displaystyle\max_{\varphi\in\Delta}}\prod_{t=1}^{m}{f(y_{t}\big{|}y_{t-1};\varphi)},$
(3)

where $\hat{\varphi}$ is the maximum likelihood estimator of ${\varphi}$.
Since the natural logarithm function increases on $(0,\infty)$, the maximum
value of the likelihood function, if it exists, occurs at the same points as
the maximum value of the logarithm of the likelihood function. In this paper,
we propose to work with the log-likelihood function which is defined as:

$\hat{\varphi}_{MLE}=\smash{\displaystyle\max_{\varphi\in\Delta}}lnL(\varphi\big{|}Y)=\smash{\displaystyle\max_{\varphi\in\Delta}}lnL_{Y}(\varphi)=\smash{\displaystyle\max_{\varphi\in\Delta}}\sum_{t=1}^{m}ln{f(y_{t}\big{|}y_{t-1};\varphi)}.$
(4)

Since this is a highly non-linear and complicated function of the unknown
parameters, we first consider the initial state vector ${\bf x}_{0}$ and
develop a set of recursions for the log-likelihood function with its first two
derivatives [12]. We then use Newton-Raphson algorithm [14] successively until
the negative of the log-likelihood is minimized to obtain the MLE.

###  2.3 GARCH model

The Generalized Autoregressive Conditional Heteroscedasticity (GARCH) model
[2, 3], was introduced in order to model the fluctuations of the variances of
financial data. It is conditional, because in these models, the nature of
subsequent volatility is conditioned by the information of the current period.
Heteroscedasticity refers to non-constant volatility. The observations $y_{t}$
of high frequency financial time series used in this paper may be represented
as:

$\displaystyle y_{t}$ $\displaystyle=\sigma_{t}\eta_{t},$ (5)

where $\sigma_{t}$ is the volatility of the observations and
$\\{\eta_{t}\\}_{t\in\mathbb{N}}$ is a Gaussian white noise sequence,
independent of $\\{\sigma_{t}\\}_{t\in\mathbb{N}}$ and
$\\{y_{t}\\}_{t\in\mathbb{N}}$. This equation can be interpreted as the
observation equation of a state space model (see subsection 2.1), whereby the
state equation is a recursive formula for the state $\sigma_{t}$:

$\displaystyle\sigma^{2}_{t}$
$\displaystyle=a_{0}+a_{1}y^{2}_{t-1}+b_{1}\sigma^{2}_{t-1},$ (6)

where $a_{0},a_{1},b_{1}\geq 0$, so that $\sigma^{2}_{t}>0$ for any values of
$y_{t}$. Eqs. (5) and (6) admit a non-Gaussian ARMA (1,1) model [13] for the
squared process as:

$y^{2}_{t}=a_{0}+(a_{1}+b_{1})y^{2}_{t-1}+\phi_{t}-b_{1}\phi_{t-1},$ (7)

where $\phi_{t}=\sigma_{t}^{2}(\eta_{t}^{2}-1)$. In order to compute the
variance at time t, we follow the standard GARCH $(m,n)$ model which is of the
form:

$\sigma^{2}_{t}=a_{0}+\sum_{j=1}^{m}a_{j}y_{t-j}^{2}+\sum_{j=1}^{n}b_{j}\sigma_{t-j}^{2}.$
(8)

If $n=0$ then the GARCH model changes into an ARCH (m) model.

The parameters $a_{0},a_{i},$ and $b_{j}$ are estimated by MLE (subsection
2.2) using the lilkelihood function. Taking into account the Normal
probability density function, the conditional likelihood in Eq. (4) is
obtained from the product of Normal ($N(0,\sigma_{t}^{2})$) densities with
$\sigma_{t}^{2}$. Using the estimated parameters, we obtain one-step-ahead
prediction of the volatility ($\widehat{\sigma}^{2}_{t}$), that is,

$\widehat{\sigma}^{2}_{t}=\widehat{a}_{0}+\sum_{j=1}^{m}\widehat{a}_{j}y_{t+1-j}^{2}+\sum_{j=1}^{n}\widehat{b_{j}}\widehat{\sigma}_{t+1-j}^{2}.$
(9)

We can analyze the residuals and squared residuals to test the Normality using
some statistical tests, for instance, Jarqua-Bera test [15], Shapiro-Wilk test
[16], Ljung-Box test [17], and LM-Arch test [18].

###  2.4 Stochastic Volatility model

The stochastic volatility (SV) model implies that the volatility is driven by
an innovation sequence, that is, independent of the observations [19]. It
causes volatility through an unobservable process that allows volatilities to
vary stochastically. To develop the SV model, we use the log-squared
observations of the time series in Eq. (5):

$logy^{2}_{t}=log\sigma^{2}_{t}+log\eta^{2}_{t}$ $\Rightarrow
g_{t}=v_{t}+log\eta^{2}_{t},$ (10)

where $g_{t}=logy^{2}_{t}$ and $v_{t}=log\sigma^{2}_{t}$. This equation is
considered as the observation equation, and the stochastic variance $v_{t}$ is
known to be an unobserved state process. Considering the autoregression, the
form of $v_{t}$ can be expressed as:

$v_{t}=\alpha_{0}+\alpha_{1}v_{t-1}+\omega_{t},$ (11)

where $\omega_{t}$ is a white Gaussian noise with variance
$\sigma_{\omega}^{2}$. Eqs. (10) and (11) constitute the stochastic volatility
model by Taylor [20]. In this study, our approach is to estimate the
parameters $\alpha_{0},\alpha_{1},\sigma_{\omega}$ and then forecast the
future observations $y_{n+m}$ from $n$ data points.

To compute the observation noise, we consider the mixtures of two Normal
distributions with one centered at zero. Thus we have:

$g_{t}=\lambda+v_{t}+\gamma_{t},$ (12)

where $\lambda$ is the mean of log-squared observations and
$\gamma_{t}=Q_{t}z_{t0}-(Q_{t}-1)z_{t1}$, which fulfills the following
condition:

$\displaystyle z_{t0}\sim\text{i.i.d}\mbox{ N}(0,\sigma^{2}_{0})$
$\displaystyle z_{t1}\sim\text{i.i.d}\mbox{ N}(\phi_{1},\sigma^{2}_{1})$
$\displaystyle Q_{t}\sim\text{i.i.d}\mbox{ Bernoulli }(p),$

where $p$ is an unknown mixing probability. We therefore define the time-
varying probabilities $\text{Pr}\\{Q_{t}=0\\}=p_{0}$ and
$\text{Pr}\\{Q_{t}=1\\}=p_{1}$, where $p_{0}+p_{1}=1$.

####  2.4.1 Parameter Estimation

In order to estimate the parameters, we use the filtering technique that is
followed by three steps namely, forecasting, updating, and parameter
estimation. In the first step, we forecast the unobserved state vector $s_{t}$
on time series observations as follows:

$v^{t}_{t+1}=\alpha_{0}+\alpha_{1}v^{t-1}_{t}+\sum_{j=0}^{1}p_{tj}K_{tj}\eta_{tj},$
(13)

where the predicted state estimators
$v^{t-1}_{t}=E(v_{t}|y_{1},\ldots,y_{t-1})$. The corresponding error
covariance matrix is defined as:

$U^{t}_{t+1}=\alpha^{2}_{1}U^{t-1}_{t}+\sigma_{\omega}^{2}-\sum_{j=0}^{1}p_{tj}K^{2}_{tj}\sum\nolimits_{tj}.$
(14)

At this point, the innovation covariances are given as
$\sum_{t0}=U_{t}^{t-1}+\sigma_{0}^{2}$ and
$\sum_{t1}=U_{t}^{t-1}+\sigma_{1}^{2}$, where
$U^{t-1}_{t}=TU^{t-1}_{t-1}T^{t}+\sigma_{\omega}^{2}$, $U^{0}_{0}=\sum_{0}$,
and $\sum_{t}=\text{var}(\eta_{t})$. Furthermore, we use Kalman filter [21] to
measure the estimates precision, which may be shown as:

$K_{t0}=\alpha_{0}U_{t}^{t-1}/(U_{t}^{t-1}+\sigma_{0}^{2})\mbox{ and
}K_{t1}=\alpha_{1}U_{t}^{t-1}/(U_{t}^{t-1}+\sigma_{1}^{2}).$ (15)

The second step deals with updating results while we have a new observation of
$y_{t}$ at time t. The prediction errors of the likelihood function are
computed using the following relations:

$\eta_{t0}=g_{t}-\lambda-v_{t}^{t-1}\mbox{ and }\eta_{t1}=g_{t}-\lambda-
v_{t}^{t-1}-\phi_{1}.$ (16)

For estimating the parameters, we complete the updating step by assessing the
time-varying probabilities (for $t=1,\ldots,m$):

$\displaystyle p_{t1}$
$\displaystyle=\frac{p_{1}h_{1}(t|t-1)}{p_{0}h_{0}(t|t-1)+p_{1}h_{1}(t|t-1)}$
$\displaystyle\text{and }p_{t0}$ $\displaystyle=1-p_{t1},$

where $h_{j}(t|t-1)$ is considered to be the conditional density of $y_{t}$,
given the previous observations $y_{1},\ldots,y_{t-1}$.

Since the observation noise of this model is not fully Gaussian, it is
computationally difficult to obtain the exact values of $h_{j}(t|t-1)$. Hence,
we use a good approximation of $h_{j}(t|t-1)$ that provides Normal density
which is: $N(v_{t}^{t-1}+\phi_{j},\sum_{tj}),$ for $j=0,1$ and $\phi_{0}=0.$

Finally, we estimate the parameters
($\Theta=(\alpha_{0},\alpha_{1},\sigma_{w}$,
$\sigma_{0},\phi_{1},\sigma_{1})^{\prime})$ by maximizing the expected
likelihood, where the MLE is represented as:

$lnL_{Y}(\Theta)=\sum_{t=1}^{m}ln\big{(}\sum_{j=0}^{1}p_{j}h_{j}(t|t-1)\big{)}.$
(17)

##  3 Dynamic behavior of the data sets

In this section, we present the background of the time series arising in
finance and geophysics. It is the dynamic behavior of the data that encourages
us to apply our methodology in this paper.

###  3.1 Financial Time Series

We study a set of high-frequency financial returns (per minute) from the
following four stock exchanges: Bank of America Corporation (BAC), Discover
Financial Services (DISCOVER), INTEL semiconductor manufacturing company
(INTEL), and IAMGOLD corporation (IAG). Fig. 1-4 provide a good perspective on
the trending direction or risk management of the high frequency returns of
stock markets. The financial crisis that occurred is evident in the large
spikes in the figures. We see that the volatility of data changes dramatically
at a short interval and that the periods of high volatility are sometimes
correlated. This establishes that the volatility itself is very volatile. The
fluctuations of stock returns (per minute) typically exhibit the volatility
clustering. That is to say, small changes in the price tend to be followed by
small changes, and large changes by large ones. The volatility clustering
suggests that current information is highly correlated with past information
at different levels.

Figure 1: Financial returns of high-frequency trading observations (per
minute) from Bank of America stock exchange. Figure 2: Financial returns of
high-frequency trading observations (per minute) from DISCOVER Financial
Services stock exchange. Figure 3: Financial returns of high-frequency trading
observations (per minute) from INTEL Corporation stock exchange. Figure 4:
Financial returns of high-frequency trading observations (per minute) from IAG
stock exchange.

###  3.2 Geophysical Time Series

The time series used in this study correspond to a set of magnitude 3.0-3.3
aftershocks of a recent magnitude 5.2 intraplate earthquake which occurred on
June 26, 2014. These earthquakes were located near the town of Clifton,
Arizona, where a large surface copper mine previously triggered off several
explosions forming part of quarry blasts activities. We selected some
explosions cataloged with similar magnitude as the earthquakes (M=3.0-3.3) and
located in the same region within a radius of 10km. We collected the
seismograms containing the seismic waves from two nearby seismic stations
(IU.TUC and IU.ANMO) located between 150 and 400km from the seismic events.
The data contains information about the date, time, longitude, latitude, the
average distance to seismic events, average azimuth, and the magnitude of each
seismic event in the region (see Tables 3.2 and 3.2). The dynamic behavior of
the geophysical time series is shown in Figs. 6-9. In these figures, we notice
that the frequency components change from one interval to another in
earthquake or explosion as long as it lasts. The mean of the series appears to
be stable with an average magnitude of approximately zero. This dynamic
behavior illustrates the time evolution of the magnitude with its volatility.
The volatility depends on time, in that it is high at a certain point, but low
at another. The volatility clustering reflects its varying nature in time, as
well as the mean reversion characteristics of the data.

Table 1: Stations information

Station | Network | Latitude | Longitude |  Average distance to events (km) |  Average Azimuth (deg)
---|---|---|---|---|---
TUC | IU | $32.3^{\circ}$ | $-110.8^{\circ}$ |  161 |  76
ANMO | IU | $34.9^{\circ}$ | $-106.5^{\circ}$ |  357 |  224

Table 2: Events information Event | Magnitude | Date | Time (UTC) | Latitude | Longitude
---|---|---|---|---|---
Earthquake | 3.0 | 7/12/14 | 7:12:53 | $32.58^{\circ}$ | $-109.08^{\circ}$
Explosion | 3.2 | 12/23/99 | 21:15:48 | $32.65^{\circ}$ | $-109.08^{\circ}$
Figure 5: The map shows the location of the seismic stations IU.TUC and
IU.ANMO used in this study (yellow color triangles). Red open circle
represents the area within which the earthquakes and explosions used in this
study are located. Figure 6: The arrival phases from an earthquake in Table-1
as recorded by TUC station. Figure 7: The arrival phases from an explosion in
Table-1 as recorded by TUC station. Figure 8: The arrival phases from an
earthquake in Table-1 as recorded by ANMO station. Figure 9: The arrival
phases from an explosion in Table-1 as recorded by ANMO station.

##  4 Stationary Approach

In this section, we analyze the time series by testing for stationarity in the
high frequency financial data and seismic waves generated by the earthquake
and explosion data by using the Augmented Dickey Fuller (ADF) [22] and
Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests. These two tests are very
powerful and capable of handling very complex models.

###  4.1 ADF test

The ADF test tests the null hypothesis that a time series $y_{t}$ is a unit
root against the alternative that it is stationary, assuming that the dynamics
in the data have an ARMA structure [23]. The summary statistics of this test
for the data sets used in this work are displayed in Tables 3 and 4
respectively.

Table 3: ADF t-statistics test for financial data Stocks | Test statistics | P-value
---|---|---
Bank of America | -14.873 | 0.01
Discover | -13.955 | 0.01
Intel | -21.314 | 0.01
IAG | -13.830 | 0.01
Table 4: ADF t-statistics test for geophysical data Events | TUC station | ANMO station
---|---|---
Test statistics | p-value | Test statistics | p-value
Earthquake | -42.018 | 0.01 | -40.509 | 0.01
Explosion | -40.831 | 0.01 | -38.954 | 0.01

Test interpretation:
$H_{0}:$ There is a unit root for the time series.
$H_{a}:$ There is no unit root for the time series. This series is stationary.

The t-statistics are used to compute the p-values, which are compared with the
significance level (0.05) and which suggest whether the null hypothesis is
acceptable or not. Since the computed p-value is lower than the significance
level $\alpha=0.05$, we reject the null hypothesis $H_{0}$ for both financial
and geophysical data, and accept the alternative hypothesis $H_{a}$. Thus the
data under study are stationary in time.

###  4.2 KPSS test

The KPSS test are used for testing a null hypothesis that an observable time
series is stationary against the alternative of a unit root [24]. The summary
statistics for the results of this test are displayed in Tables 5 and 6
respectively.

Table 5: KPSS t-statistics test for financial data Stocks | Test statistics | P-value
---|---|---
Bank of America | 0.1733 | 0.1
Discover | 0.1221 | 0.1
Intel | 0.0069 | 0.1
IAG | 0.3412 | 0.1
Table 6: KPSS t-statistics test for geophysical data Events | TUC station | ANMO station
---|---|---
Test statistics | p-value | Test statistics | p-value
Earthquake | 0.0017 | 0.1 | 0.0020 | 0.1
Explosion | 0.0012 | 0.1 | 0.0030 | 0.1

Test interpretation:
$H_{0}:$ There is a unit root for the time series. This series is stationary.
$H_{a}:$ There is no unit root for the time series.

As the computed p-value is greater than the significance level $\alpha=0.05$,
we accept the null hypothesis $H_{0}$ in all data sets. Thus the time series
used in this paper are stationary time series.

In the next two sections we describe the analysis of time series arising in
geophysics and finance using a deterministic and stochastic approach. First we
present the estimates when the models are applied to the data sets. To
estimate the time-varying parameters, we used a deterministic and stochastic
model on the time series. The idea is to compare the two techniques and
observe which is suitable for forecasting the volatility. The analysis was
performed by an R statistical software module.

##  5 Analysis of the Deterministic Model

We present the results of the estimated parameters of high frequency financial
data obtained with the GARCH model. The red line in the following figures show
the theoretical probability density function of Normal distribution with the
same mean and standard deviation as the financial data. We therefore consider
the ARCH Normality assumption on the basis of volatility $\eta_{t}$.

(a) BAC (b) DISCOVER (c) INTEL (d) IAG

Figure 10: The histograms of financial time series and the fitted Normal
density.

The estimates of the parameters $a_{0},a_{1}$, and $b_{1}$ of the GARCH model
are stable, as the GARCH-statistic shows ( see Tables 7-8). Also, the
estimated standard errors of the parameters in most cases are small. The
smaller p-values ($<$ significance level) provide strong evidence that the
GARCH $(1,1)$ model with the specified parameters is a good fit for our data.
The volatility level of persistence can be determined by non-negative
parameters $a_{1}$ and $b_{1}$ from these tables. We can see that the
constraint $(a_{1}+b_{1})$ is less than 1, which is consistent with the
existence of a stationary solution, and supports the results of stationary
tests in section 4.

Tables 11-14 summarize the standardized residuals (R) tests for Bank of
America, Discover Financial Services, INTEL semiconductor manufacturing
company, and IAG stock exchanges respectively. The Jarque-Bera and Shapiro-
Wilk tests of Normality strongly reject the null hypothesis that the white
noise innovation process ${\eta_{t}}$ is Gaussian. The p-values $(>0.05)$ of
Ljung-Box test for squared residuals (at lag $10,15,20$) and LM-Arch test
suggest that the model fits the data well, with the exception of the non-
normality of ${\eta_{t}}$. It is because the null hypothesis cannot be
rejected at any reasonable level of significance.

To facilitate the understanding of forecasting concepts, we superimpose the
plot of one-step-ahead predicted volatility and $\pm 2$ standard prediction
errors in Figs. 11-14. The predicted volatility with $\pm
2\widehat{\sigma_{t}}$ is displayed as a dashed line surrounding the original
output. It visually shows how values of predicted volatility differ over time.

We also see some limitations of the GARCH (1,1) model itself, when applied to
our financial data sets. The positive and negative high frequency returns have
the same effect, because volatility depends on squared returns. Thus it does
not help to understand the source of variations of a financial time series,
i.e. the causes of the variation in volatility. The model provides only a
mechanical way to describe the behavior of conditional variance. The blue
dashed lines during the financial crash (in Figs. 11-14) indicate that the
model tends to arbitrarily predict volatility, because it slowly responds to
large isolated returns.

Table 7: GARCH statistics for BAC stock exchange Parameter | Estimate | Error | t-statistic | p-value
---|---|---|---|---
$a_{0}$ | 1.25E-07 | 1.42E-08 | 8.788 |  $<$2E-06
$a_{1}$ | 2.47E-01 | 5.91E-02 | 4.171 | 3.0E-05
$b_{1}$ | 7.21E-01 | 3.63E-02 | 19.84 |  $<$2E-16
Table 8: GARCH statistics for the DISCOVER stock exchange Parameter | Estimate | Error | t-value | p-value
---|---|---|---|---
$a_{0}$ | 2.86E-08 | 3.89E-09 | 7.349 | 2.0E-13
$a_{1}$ | 1.79E-01 | 2.42E-02 | 7.401 | 1.4E-13
$b_{1}$ | 7.09E-01 | 3.14E-02 | 22.58 |  $<$2E-16
Table 9: GARCH statistics for INTEL stock exchange Parameter | Estimate | Error | t-statistic | p-value
---|---|---|---|---
$a_{0}$ | 0.02167 | 0.00982 | 2.206 | 0.0274
$a_{1}$ | 0.15914 | 0.03682 | 4.322 | 1.5E-05
$b_{1}$ | 0.61185 | 0.13686 | 4.471 |  $<$7.8E-06
Table 10: GARCH statistics for IAG stock exchange Parameter | Estimate | Error | t-statistic | p-value
---|---|---|---|---
$a_{0}$ | 1.31E-07 | 1.59E-08 | 8.278 | 2.2E-16
$a_{1}$ | 2.77E-01 | 4.03E-02 | 6.867 | 6.6E-12
$b_{1}$ | 4.83E-01 | 5.24E-02 | 9.208 |  $<$2E-16
Table 11: Standardised Residuals Tests for BAC stock exchange | Residuals | Tests | Statistics | p-value
---|---|---|---|---
Jarque-Bera Test | $R$ | $\chi^{2}$ | 4918340 | 0
Shapiro-Wilk Test | $R$ | W | 0.649712 | 0
Ljung-Box Test | $R$ | Q(10) | 14.71311 | 0.1429
Ljung-Box Test | $R$ | Q(15) | 17.82156 | 0.2722
Ljung-Box Test | $R$ | Q(20) | 18.77450 | 0.5365
Ljung-Box Test | $R^{2}$ | Q(10) | 0.107465 | 1
Ljung-Box Test | $R^{2}$ | Q(15) | 0.130187 | 1
Ljung-Box Test | $R^{2}$ | Q(20) | 0.147051 | 1
LM-Arch Test | $R$ | T$R^{2}$ | 0.126933 | 1
Table 12: Standardised Residuals Tests for DISCOVER stock exchange |  |  | Statistics | P-value
---|---|---|---|---
Jarque-Bera Test | R | $\chi^{2}$ | 1688779 | 0
Shapiro-Wilk Test | R | W | 0.722888 | 0
Ljung-Box Test | R | Q(10) | 10.93070 | 0.3629
Ljung-Box Test | R | Q(15) | 13.58973 | 0.5568
Ljung-Box Test | R | Q(20) | 16.81071 | 0.6652
Ljung-Box Test | $R^{2}$ | Q(10) | 0.198119 | 0.9999
Ljung-Box Test | $R^{2}$ | Q(15) | 0.258743 | 1
Ljung-Box Test | $R^{2}$ | Q(20) | 0.369189 | 1
LM Arch Test | R | $TR^{2}$ | 0.208586 | 1
Table 13: Standardised Residuals Tests for INTEL stock exchange | Residuals | Tests | Statistics | p-value
---|---|---|---|---
Jarque-Bera Test | $R$ | $\chi^{2}$ | 1769.854 | 0
Shapiro-Wilk Test | $R$ | W | 0.971157 | 0
Ljung-Box Test | $R$ | Q(10) | 296.0897 | 0
Ljung-Box Test | $R$ | Q(15) | 301.6321 | 0
Ljung-Box Test | $R$ | Q(20) | 304.7640 | 0
Ljung-Box Test | $R^{2}$ | Q(10) | 7.968585 | 0.6319
Ljung-Box Test | $R^{2}$ | Q(15) | 16.48274 | 0.3507
Ljung-Box Test | $R^{2}$ | Q(20) | 59.49665 | 8.5E-06
LM-Arch Test | R | T$R^{2}$ | 14.76281 | 0.2547
Table 14: Standardised Residuals Tests for IAG stock exchange | Residuals | Tests | Statistics | p-value
---|---|---|---|---
Jarque-Bera Test | $R$ | $\chi^{2}$ | 126041.2 | 0
Shapiro-Wilk Test | $R$ | W | 0.842114 | 0
Ljung-Box Test | $R$ | Q(10) | 12.66334 | 0.2431
Ljung-Box Test | $R$ | Q(15) | 15.02200 | 0.4498
Ljung-Box Test | $R$ | Q(20) | 21.73719 | 0.3549
Ljung-Box Test | $R^{2}$ | Q(10) | 0.146992 | 1
Ljung-Box Test | $R^{2}$ | Q(15) | 0.333308 | 1
Ljung-Box Test | $R^{2}$ | Q(20) | 0.706543 | 1
LM-Arch Test | R | T$R^{2}$ | 0.243845 | 1

GARCH predicted volatility

Figure 11: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from BAC stock exchange.

GARCH predicted volatility

Figure 12: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from DISCOVER stock exchange.

GARCH predicted volatility

Figure 13: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for one hundred observations from INTEL stock exchange.

GARCH predicted volatility

Figure 14: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from IAG stock exchange.

##  6 Analysis of Stochastic Model

In this section we use the stochastic volatility model to forecast the
volatility of a time series. The first subsection is devoted to analyzing high
frequency financial time series and the latter subsection includes geophysical
time series.

The GARCH model used in the previous section differs from the SV model because
it does not have any stochastic noise. The SV model is characterized by the
fact that it invariably contains its probability density function. We compute
the maximum likelihood by taking into consideration the conditional Normal
distribution. From Figs. 10 and 15 we conclude that the histograms of the
financial and geophysical time series respectively, are well represented.

(a) Earthquake (TUC station) (b) Explosion (TUC station) (c) Earthquake (ANMO
station) (d) Explosion (ANMO station)

Figure 15: The histograms of geophysical time series and the fitted Normal
density.

The parameters from time-varying Eqs. (11) and (12) were initialized in order
to observe the performance of the SV algorithms during a set of magnitudes for
each seismic event. We set the initial values to be
$\alpha_{0}=0,\alpha_{1}=0.96,\sigma_{\omega}=0.3,\sigma_{0}=1,\phi_{1}=-4,\sigma_{1}=3$
and $\lambda$, the mean of the observations. In order to maximize Eq. (17),
the innovation processes for Eqs. (11) and (12) were fitted to the data by
taking into consideration this time-varying probability ($p_{1}=0.5$).

###  6.1 Volatility of Financial Time Series

Using the stochastic volatility model, we estimate the parameters and one-
step-ahead predicted log-volatility of high frequency returns from four stock
exchanges.

Table 15: Summary statistics for DISCOVER stock exchange Parameter | Estimate | Standard Error
---|---|---
$\alpha_{0}$ | 0.0053 | 0.0008
$\alpha_{1}$ | 0.9799 | 0.0036
$\sigma_{\omega}$ | 0.1841 | 0.0004
$\lambda$ | -16.801 | 0.1687
$\sigma_{0}$ | 1.0697 | 0.0008
$\phi_{1}$ | -2.2973 | 0.0011
$\sigma_{1}$ | 2.8771 | 0.0008
Table 16: Summary statistics for INTEL stock exchange Parameter | Estimate | Standard Error
---|---|---
$\alpha_{0}$ | -0.1650 | 0.4123
$\alpha_{1}$ | 0.8629 | 0.0682
$\sigma_{\omega}$ | 0.5978 | 0.1871
$\lambda$ | -14.166 | 2.9146
$\sigma_{0}$ | 0.9478 | 0.1311
$\phi_{1}$ | -3.3405 | 0.1904
$\sigma_{1}$ | 3.2604 | 0.1158
Table 17: Summary statistics for IAG stock exchange Parameter | Estimate | Standard Error
---|---|---
$\alpha_{0}$ | -0.1994 | 0.4848
$\alpha_{1}$ | 0.8440 | 0.0738
$\sigma_{\omega}$ | 0.6462 | 0.1976
$\lambda$ | -14.164 | 3.0113
$\sigma_{0}$ | 0.9536 | 0.1434
$\phi_{1}$ | -3.3967 | 0.2038
$\sigma_{1}$ | 3.6077 | 0.1261
Table 18: Summary statistics for BAC stock exchange Parameter | Estimate | Standard Error
---|---|---
$\alpha_{0}$ | -0.1596 | 0.4066
$\alpha_{1}$ | 0.8641 | 0.0676
$\sigma_{\omega}$ | 0.5949 | 0.1856
$\lambda$ | -14.191 | 2.9078
$\sigma_{0}$ | 0.9471 | 0.1299
$\phi_{1}$ | -3.3374 | 0.1895
$\sigma_{1}$ | 3.238 | 0.1152

Predicted log-volatility of DISCOVER stock exchange

Figure 16: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for DISCOVER stock exchange.

Predicted log-volatility of INTEL stock exchange

Figure 17: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for INTEL stock exchange.

Predicted log-volatility of IAG stock exchange

Figure 18: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for IAG stock exchange.

Predicted log-volatility of BAC stock exchange

Figure 19: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for BAC stock exchange.

###  6.2 Volatility of Geophysical Time Series

Similarly, using the stochastic volatility model we estimate the parameters
and one-step-ahead predicted log-volatility of earthquake and explosion time
series.

Table 19: Summary statistics for Earthquake data sets Parameter | TUC station | ANMO station |
---|---|---|---
Estimates | Standard Error | Estimates | Standard Error |
$\alpha_{0}$ | 0.0286 | 0.0343 | 0.0828 | 0.1015 |
$\alpha_{1}$ | 0.9960 | 0.0013 | 0.9851 | 0.0028 |
$\sigma_{\omega}$ | 0.4534 | 0.0357 | 0.7427 | 0.0183 |
$\lambda$ | -4.0974 | 7.4020 | -5.2073 | 6.4241 |
$\sigma_{0}$ | 0.6284 | 0.0401 | 0.0007 | 0.0720 |
$\phi_{1}$ | -2.4730 | 0.0862 | -2.3090 | 0.0773 |
$\sigma_{1}$ | 2.3481 | 0.0518 | 2.1530 | 0.0482 |
Table 20: Summary statistics for Explosion data sets Parameter | TUC station | ANMO station |
---|---|---|---
Estimates | Standard Error | Estimates | Standard Error |
$\alpha_{0}$ | 0.1778 | 0.1025 | 0.1889 | 0.1121 |
$\alpha_{1}$ | 0.9874 | 0.0028 | 0.9814 | 0.0032 |
$\sigma_{\omega}$ | 0.7492 | 0.0181 | 0.7003 | 0.0171 |
$\lambda$ | -10.151 | 6.7569 | -10.076 | 5.2518 |
$\sigma_{0}$ | 2.13E-05 | 0.0517 | 7.09E-05 | 0.0501 |
$\phi_{1}$ | -2.3475 | 0.0798 | -2.3063 | 0.0740 |
$\sigma_{1}$ | 2.1657 | 0.0491 | 2.0917 | 0.0462 |

Predicted log-volatility of Earthquake data

Figure 20: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from TUC station.

Predicted log-volatility of Explosion data

Figure 21: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from TUC station.

Predicted log-volatility of Earthquake data

Figure 22: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from ANMO station.

Predicted log-volatility of Explosion data

Figure 23: One-step-ahead predicted log-volatility, with $\pm 2$ standard
prediction errors for two hundred observations from ANMO station.

###  6.3 Discussion

Tables 17 to 20 summarize the estimation of parameters
($\alpha_{0},\alpha_{1}$, $\sigma_{w}$, $\lambda$, $\sigma_{0}$, $\phi_{1}$,
and $\sigma_{1}$). The estimated error in these tables makes two things
evident: firstly, the estimates are close to the true parameters; secondly,
the algorithm of the SV model is consistent with the results obtained by using
the data. The variance $\sigma_{w}^{2}$ of the log-volatility process measures
the uncertainty related to the future data volatility. If the value of
$\sigma_{w}^{2}$ is zero, it is not possible to identify the SV model. The
parameter $\alpha_{1}$ is considered as a measure of the persistence of shocks
to the volatility. Tables 17 to 20 indicate that $\alpha_{1}$ is less than 1,
which suggests that the latent volatility process and $y_{t}$ are stationary,
what confirms the results of section 4.

In these tables, we notice that $\alpha_{1}$ is near to unity and
$\sigma_{w}^{2}$ is different from 0, which means that the volatility
evolution is not smooth over time. It also suggests that the time series could
be heteroscedastic by nature, that is, there is a non-constant conditional
volatility over time. So, it is very useful to control the risk or to mitigate
the effect of hazards. For example, if there are two time series having the
same mean but with different variances, we would then consider the series with
lower variance, because it is less risky.

##  7 Conclusion

In this study, we have used high frequency financial returns and geophysical
time series measured every minute. We have implemented a type of volatility
models that incorporates time-varying parameters in a stationary scenario. To
obtain a good fit for the data, we used a deterministic model for high-
frequency financial returns and a stochastic model for the financial and
geophysical time series. The geophysical data aligns with the SV model because
of its stochastic behavior whereas the GARCH model does not fit the
geophysical data. This is because the GARCH model tends to arbitrarily predict
volatility in the case of large arrival phases of earthquake and explosion
time series.

The parameter estimates of the GARCH (1,1) model indicate that there exists a
stationary solution in the conditional volatility of high frequency financial
returns (see Section 5). For the SV model, we estimated the volatility
parameters of time series based on the geophysical and financial time series.
The fitted SV model allows us to capture the volatility evolution that
suggests the physical and long-memory behavior of the data. With the use of
maximum likelihood computation, we have succeeded in making a good prediction
despite the variation of the observational noise from a Normal mixture
distribution, because the geophysical time series studied is not fully
Gaussian (see the histograms in Fig. 15).

Our results suggest that the stochastic process to forecast the time series is
more effective in enforcing the characteristic of time-varying parameters than
the commonly used deterministic process. It is because the one-step-ahead
predictions along with the estimated standard error of stochastic volatility
model do not show any limitations unlike the GARCH model (see section 5). We
notice that the GARCH model is more sensitive to noise or unexpected shocks
than the stochastic volatility model. The advantage of the stochastic
methodology is that the estimates obtained are stable around the true value.
Moreover, the low errors (see Tables 19 and 20) imply that the estimation
procedure is accurate, meaning that it is capable of generating a higher
forecasting accuracy.

## References

  * [1] Fong, S.J., and Nannan, Z. (2011), Towards an Adaptive Forecasting of Earthquake Time Series from Decomposable and Salient Characteristics, _The Third International Conferences on Pervasive Patterns and Applications - ISBN : 978-1-61208-158-8_ , 53-60.
  * [2] Engle, R.F. (1982), Autoregressive Conditional Heteroskedasticity with Estimates of the Variance of United Kingdom Inflation, _Econometrica_ , 50(4), 987-1007.
  * [3] Bollerslev, T. (1986), Generalized Autoregressive Conditional Heteroskedasticity, _J. Econometrics_ , 31, 307-327.
  * [4] Hamiel, Y., Amit, R., Begin, Z.B., Marco, S., Katz, O., Salamon, A., Zilberman, E., and Porat, N. (2009), The seismicity along the Dead Sea fault during the last 60,000 years. _Bulletin of Seismological Society of America_ ,  99(3), 2020-2026.
  * [5] Mariani, M.C. and Tweneboah, O.K. (2016), Stochastic differential equations applied to the study of geophysical and financial time series, _Physica A_ , 443, 170-178.
  * [6] Mariani, M.C., Bhuiyan, M.A.M, and Tweneboah, O.K. (2018), Estimating Stochastic Volatility by using Ornestein-Uhlenbeck Type Models, _Physica A_ \- Statistical Mechanics and its Applications, 491, 167-176.
  * [7] Brockman, P., and Chowdhury, M. (1997), Deterministic versus stochastic volatility: implications for option pricing models, _Applied Financial Economics_ , 7, 499-505.
  * [8] Mariani, M.C., Gonzalez-Huizar, H., Bhuiyan, M.A.M. and Tweneboah, O.K.T. (2017), Using Dynamic Fourier Analysis to Discriminate Between Seismic Signals from Natural Earthquakes and Mining Explosions, _AIMS Geosciences_ , 3(3), 438-449.
  * [9] Rubio, F. J., and Johansen, A. M. (2013), A simple approach to maximum intractable likelihood estimation, _Electronic Journal of Statistics_ , 7, 1632-1654.
  * [10] Commandeur, J.F., and Koopman, S.J. (2007), An Introduction to State Space Time Series Analysis, _Oxford University press_ , 107-121.
  * [11] Eliason, S.R. (1993), Maximum Likelihood Estimation-Logic and Practice , _Quantitative applications in the social sciences_ , 96, 1-10.
  * [12] N. K. Gupta and R. K. Mehra (1974), Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations, _IEEE Transactions on Automatic Control_ , 19(6), 774-783.
  * [13] Rajarshi, M.B. (2012), Statistical Inference for Discrete Time Stochastic Processes, _Springer - ISBN: 978-81-322-0763-4_ , 39-55.
  * [14] Jones, R.H. (1980), Maximum likelihood fitting of ARMA models to time series with missing observations, _Technometrics_ , 22(3), 389-395.
  * [15] Brys, G., Hubert, M., and Struyf, A. (2004), A Robustification of the Jarque-bera test of normality, _Physica-Verlag/Springer_ , 753-760.
  * [16] Shapiro, S.S., and Wilk, M.B. (1965), An Analysis of Variance Test for Normality (Complete Samples), _Biometrika_ , 52(3/4), 591-611.
  * [17] Ljung, G.M., and Box, G.E.P. (1978), On a Measure of Lack of Fit in Time Series Models, _Biometrika_ , 65(2), 297-303.
  * [18] Wang, W., Gelder, P.H.A.J.M.V., Virijling, J.K., and Ma, J. (2005), Testing and modelling autoregressive conditional heteroskedasticity of streamflow processes, _European Geosciences Union_ , 12(1), 55-66.
  * [19] Janssen, A., and Drees, H. (2016), A stochastic volatility model with flexible extremal dependence structure, _Bernoulli_ , 22(3), 1448-1490.
  * [20] Taylor, S. J. (1982). Financial returns modeled by the product of two stochastic processes, A study of daily sugar prices, 1961-79. _Time Series Analysis: Theory and Practice, ZDB-ID 7214716_ , 1, 203-226.
  * [21] Cipra, T. and Romera, T. (1991), Robust Kalman Filter and Its Application in Time Series Analysis, _Kybernetika_ , 27(6), 481-494.
  * [22] Mariani, M.C., Huizar, Bhuiyan, M.A.M., Tweneboah, O.K. (2017)- Using Dynamic Fourier Analysis to Discriminate Between Seismic Signals from Natural Earthquakes and Mining Explosions, _AIMS Geosciences_ , 3(3), 438-449.
  * [23] Said, S.E. and Dickey, D. (1984), Testing for Unit Roots in Autoregressive Moving-Average Models with Unknown Order, _Biometrika_ , 71, 599-607.
  * [24] Kwiatkowski, D., Phillips, P.C.B., Schmidt, P., and Shin, Y. (1992), Testing the null hypothesis of stationarity against the alternative of a unit root, _Journal of Econometrics_ , 54, 159-178.