random sample ˜xt+1:T drawn approximately from p (xt+1:T|z1:T), take one step back in time and sample ˜xt from p (xt|˜xt+1, z1:T). Repeating this process sequentially back over time produces the smoothing algorithm.
6.5 Parameter Estimation using EM Algorithm with Particle Filtering Method
In this section, we take SV-DEJ-VIJ model as example, let {rt}Tt=1 be the observable log return. The jump number of the return {Ny,tP }Tt=1, the upward jump size of return {ξy,t+}Tt=1, the downward jump size of the return {ξ−y,t}Tt=1, the stochastic volatility {vt}Tt=0, jump number of the return {Nv,t}Tt=1, and the jump size of the volatility {ξv,t}Tt=1 are the latent data.
The discrete-time version of the dynamic process of SV-DEJ-VIJ model:
rt=
where ∆ is the time interval, µ is the drift term, the compound Poisson process of return jump,
1, with probability p,
−1, with probability q,
where Ny,tP is used to depict the number of the return jump which follows a Poisson process with the intensity λy, i.e., Ny,tP ∼ P oisson (λ∆), Nt is the jump direction, i.e.,
‧
compensator of the return jump and the average jump amplitude.MJPy = λy
p
1 − η+ + q
1 + η− − 1
The stochastic volatility vt follows a mean-reverting square root process with long term level θ > 0, the speed of adjustment κ > 0, and the variation coefficient of the diffusion volatilityσv > 0. The compound Poisson process of volatility jump,
Jv,t= ξv,tNv,t=
Nv,t
X
j=1
ξjv,t
where the volatility jump size ξv,t is a exponential random variable with mean µv, Nv,t is used to depict the number of the volatility jump which follows a Poisson process with the intensity λv. Moreover, εP1,t and εP2,t are independent standard normal random variables.
By the Cholesky decomposition, they are generated from the dependent Browian motions with correlation coefficient ρ. The equation (6.13) is the measure equation of the return rt which reflects the stochastic volatility and the jump risk. In order to use the particle filtering method to track the dynamic of the latent variable, we substitute the measure equation (6.13) into the equation (6.14) and set ∆ = 1, i.e., one day interval. Then we get the state (or transition) equation associating with the return rt as follows:
vt= vt−1+ κ (θ − vt−1)
where rt is the observable log return The jump number of the return Ny,tP , the upward jump size of return ξy,t+, the downward jump size of the return ξ−y,t, the stochastic volatility {vt}Tt=0, jump number of the return {Nv,t}Tt=1, and the jump size of the volatility {ξv,t}Tt=1 are the latent data. According to the state equation and measure equation, we can con-struct the state space of SV-DEJ-VIJ Model. Then we can use the Particle filtering method to track the latent variables and estimate the parameters of the SV-DEJ-VIJ model with the expectation-maximization (EM) algorithm. We divided the estimation
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
process into three steps:
Step 1. Particle Filtering Method
Given observable return rt, for t = 1, ..., T , we sample each random variables of SV-DEJ-VIJ with R particles from their corresponding distributions, that is, Ny,tP ∼ P oisson λPy, ξy,t+ ∼ e exp (η+), ξy,t− ∼ exp (η−), Nv,t ∼ P oisson (λv), ξv,t ∼ exp (µv) and εP2,t ∼ N (0, 1) . Then we can evaluate the weightcorresponding toand latent random variables through the equation (6.12) and (6.13), the weight is given by:
˜
ωt ∝ p (rt|vt) = (2πvt)−12 exp
− 1 2vt
rt−
µ − 1
2vt−1− MJPy
−
Ny,tP
X
j=1
ξy,tj
2
(6.16) Given the sampling latent random variables and the weight ˜ωt at time t, the posterior distribution of stochastic volatility is given by normalizing the weight ωt
p vit|rt = ωti = ˜ωit/
R
X
j=1
˜
ωtj, i = 1, ..., R,
Thus, we can use the resampling algorithm to resample the state variables Ny,tP, ξy,t+, ξy,t−, Nv,t and ξv,t with corresponding weight ωt.
Step 2. Smoothing
After repeating the iteration for particle filtering from time 1 to T , we obtain the particle series of the state variables. However, the particles may loss some information if a time series of the sample is a longer period of time after resampling the particles. That is to say, the resampling will cause the particles to be more monotonous and the particles can-not track the series of the state variables accurately. In order to modify those drawback, Godsill, Doucet, and West (2004) provide the backward simulation method to resample the particles in the reverse-time direction conditioning on future states. The weight for the smoothing process is associated with the equation (6.15) as follows
eωt−1|t∝ f (vt−1|vt, rt) = (2πΛ1)−12 exp
− 1
2Λ2 (vt− Λ2)2
(6.17)
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
where
Λ1 = σ2vvt−1 1 − ρ2 , Λ2 = ρσv
rt−
µ − 1
2vt−1− MJPr
+
Ny,tP
X
j=1
ξr,tj
,
Λ3 = vt−1− κ (θ − vt−1) +
Nv,tP
X
j=1
ξr,tj − Λ2
Given the sampling latent random variables and the weight ˜ωt−1|t at time t − 1, the posterior distribution of stochastic volatility is given by normalizing the weight ˜ωt−1|t
Pr(vt−1i |vit, rt) = ωit−1|t= ωet−1|ti PR
j=1eωt−1|tj , i = 1, ..., R.
Thus, we can use the smoothing algorithm in the section 6.4 to resample the state vari-ables Ny,tP , ξy,t+, ξy,t−, Nv,t and ξv,t with corresponding weight ωt−1|t from T to 1.
Step 3. EM Algorithm
After the particle filtering and the smoothing steps, we obtain the state variable of SV-DEJ-VIJ, i.e., the jump number of the return {Ny,t}Tt=1, the upward jump number of the return {Ny,t+}Tt=1, the downward jump number of the return {Ny,t−}Tt=1, the total jump size of return {Jy,t}Tt=1, the total upward jump size of return {Jy,t+}Tt=1, the total downward jump size of the return {Jy,t−}Tt=1, the stochastic volatility {vt}Tt=0, jump number of the volatility {Nv,t}Tt=1, and the total jump size of the volatility {Jv,t−}Tt=1. Now, we use the expectation-maximization (EM) algorithm to estimate parameters of SV-DEJ-VIJ model.
Let Ψ = {µ, κ, θ, σ, ρ, λy, λv, µv, p, η+, η−} be the set of the parameters of SV-DEJ-VIJ.
‧
Under the SV-DEJ-VIJ model, the complete likelihood function is L Ψ|y1:T, v0:T, Ny,1:T, Ny,1:T+ , Jy,1:T+ , Jy,1:T− , Nv,1:T, Jv,1:T
There are two steps for the EM algorithm. (i) E-step: Compute the expectation of the complete log-likelihood function given the observable data (log-return data), latent variables and the k-th estimated parameters, denoted by Q Ψ|Ψ(k).
Q Ψ|Ψ(k)
where R is the number of the particles. (ii) M-step: we find the specific parameters by maximize the equation (6.19) and obtain the k-th estimated parameters, Ψ(k+1). That is,
Ψ(k+1) = arg max
Ψ
Q Ψ|Ψ(k)
Repeating above two steps of the EM algorithm until the log-likehood function converge and obtain the estimated parameters of SV-DEJ-VIJ model. The parameters of the other models are estimated by the similar method.
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
6.6 Joint Estimation
When we calibrate the option prices using the sum of squared pricing errors (SSE) only, there is a problem with the estimated parameters which does not have significantly economic meaning. Take a parameter κ as example, under the risk neutral measure, κQ = κ + h2, the optimize algorithm only find the value κQ after calibrating. That is to say, we cannot capture the economic meaning of and risk premium h2 accurately. There-fore, in order to solve this problem, we introduce the joint estimation method using a weighted joint likelihood of weekly S&P500 index call options on each Wednesday from January 1, 2007 to August 31, 2017 and daily S&P500 index returns from January 1, 2005 to August 31, 2017 in Ornthanalai (2014). This method can improve the stability of the estimated parameters under the physical measure, and then the optimize algo-rithm can tune the risk premium exactly. However, we do some slight changes to the joint estimation method from Ornthanalai (2014). We use the rolling time-window for S&P500 index returns to do the particle filtering and smoothing and then estimate the parameters with EM algorithm. For example, on August 31, 2017, we do the particle filtering and smoothing to obtain the latent variables during the period from August 31, 2015 to August 30, 2017 and then estimate the parameter with call options in the cross section on August 31, 2017. Following Hu and Zide (2002) and Ornthanalai (2014), the weighted joint log-likelihood for each day:
LJ ointt = T1+ Nt 2
LPreturn,t
T1 +T1+ Nt 2
LQoption,t
Nt , t = 1, ..., T2, (6.20) where T1is two years, T2is the number of days in the rolling windows on each trading date, and Nt is the total number of option contracts on each trading date time t. 12(T1 + N ) serves as a scaling constant and does not impact the parameter estimates.
The log-likelihood for daily S&P500 index returns LPreturn,t is related to particle filtering, smoothing and EM algorithm at time t which is discussed in section 6.5. However, there is a issue about the log likelihood function LPreturn,t in (6.18) when we compare the model’s performance of the dynamic fitting, the log likelihood function of each model has different numbers of parameters. That is to say, we cannot directly compare the models with log likelihood function in (6.18), because they have different dimension of model’s parameters. Therefore, we change the log likelihood function which is corresponding to
‧
Note that this log likelihood function under the measure in (6.21) is only used to examine the performance of the model. When we optimize the joint likelihood function, we still use the log likelihood function in (6.18).
Moreover, in order to computer the log-likelihood for fitting call options LQoption,t, we assume that the option pricing errors related to the implied volatility follow a normal distribution,
εt,j = IVt Ot,jBS − IVt Ot,jM odel , j = 1, ..., Nt,
Assume that the implied volatility error εt,j ∼ N 0, σ2ε,t, where σε,t2 is the variance of the Black-Scholes implied volatility of the market option price for each day. Thus, the complete log-likelihood of option price:
LQoption,t = ln
where IVt Ot,jBS is the Black-Scholes implied volatility of the j-th market-observed call op-tion price Ot,jBS = C (St,j, Kt,j, τt,j, rt,j), IVt OM odelt,j is the Black-Scholes implied volatility of the j-th call option price Ot,jM odel = C (Ψt|St,j, Kt,j, τt,j, rt,j) which is computed using the model, where the parameters St,j, Kt,j, τt,j, rt,j are the underlying S&P500 index price, strike, days-to-expiration, riskless rate and Ψt is the parameter vector containing the model’s risk-neutral parameter and risk premiums ht,1 and ht,2.