Estimating parameters, the Heston model, and more…

Posted on December 8, 2010 by beggsmontecarlo

In this post, we will estimate the parameters defining a CIR process from a (simulated CIR) dataset.  We will then use antithetics to estimate an expectation dealing with the Heston model, where the underlying volatility is modeled by a CIR model.

Let X=(X_t: t \in [0,T]) be a stochastic process such that dX_t=\beta(X_t - \gamma )dt + \sigma \sqrt{X_t}dB_t, where dB_t is a standard Brownian motion.  (Note: this is the CIR process, commonly used to model interest rates.)

We arbitrarily select values for the parameters: \beta = -1, \gamma = 1, \sigma = 1, X_0=.1, T=1.

Now, we run set.seed(19870210) and generate a sample path of X_t in R, storing the simulated values for X_t in mydata.

R Code

Henceforth, we pretend that the process (X_t) with some specific but unknown parameters has generated the data.

FACT: \sum_{k}(X_{k/2^n}-X_{(k-1)/2^n})^2 \approx \int_0^T(\sigma\sqrt{X_t})^2dt

Using this, we write a function to estimate \sigma from the data.  We obtain \hat{\sigma} = 0.98

R Code

Next, let Y_t = h(X_t), where h is the Lamperati transformation for (X_t).  That is: h(x) = \int \frac{1}{\sigma\sqrt{x}}.  Then, (Y_t) is a stochastic process such that

\begin{array}{lcl} dY_t & = & \left(\frac{\beta Y_t}{2}-\frac{\beta\gamma}{\sigma^2Y_t}-\frac{1}{2Y_t}\right)+dB_t\\& = & \mu(Y_t)dt + dB_t\end{array}

where dB_t is a standard Brownian motion, and \mu is a function.  Its density (for the interval [0,T]) is given by

L_{T,Y} = \exp \left (\int_0^T\mu(Y_t)dY_t - \frac{1}{2}\int_0^T\mu^2(Y_t)dt\right),

Now define a function a(y) as the following:

\begin{array}{lcl} a(y) & = & \int \mu(y)dy \\& = & \frac{\beta y^2}{4}-\frac{\beta\gamma}{\sigma^2}\log{y}-\frac{1}{2}\log{y}\end{array}

Then:

\log(L_{T,Y} ) = a(Y_T) - a(Y_0) - \frac{1}{2}\int_0^T\mu^2(Y_t)+\mu'(Y_t)dt.

It then follows that the density L_{T,X} of (X_t) is specified by:

\log(L_{T,X} ) = a(h(X_T)) - a(h(X_0)) - \frac{1}{2}\int_0^T\mu^2(h(X_t))+\mu'(h(X_t))dt.

We can now estimate the parameters (\beta,\gamma) in R by maximizing the likelihood L_{T,X}.

R Code

Running the code, we estimate (\beta, \gamma) to be (\beta_s, \gamma_s) = (-2.42, .2).  Remember that these values are inferred from the data set we created by creating a realization of a CIR process.  Since we know the real values of the underlying process to be (\beta, \gamma) = (-1,1), we expect that if we ran this experiment n (large) times (i.e. simulate the CIR process n times and estimate (\beta, \gamma) n times), \frac{1}{n}\sum_{i} (\beta_i, \gamma_i) \approx (\beta, \gamma).

We will reference the following paper on simulating a Heston model, where the underlying volatility (instead of the interest rate is modeled by a CIR model.

A paper

In particular, consider the following basic Milstein scheme:

(10)      \hat{x}(t + \triangle) = \hat{x}(t) + \kappa(\theta - \hat{x}(t))\triangle + \epsilon\sqrt{\hat{x}(t)}Z\sqrt{\triangle} + \frac{1}{4}\epsilon^2\triangle(Z^2-1)

Claim: this scheme has a positive probability of generating negative values of \hat{x} and therefore cannot be used without suitable modifications.

Proof: Let \hat{x} = 0.  Then,

\begin{array}{lcl} \hat{x}(\triangle)& = & \kappa\theta\triangle + \frac{1}{4}\triangle(Z^2-1)\\& = & \frac{\triangle}{4}(4\kappa\theta +Z^2-1)\end{array}

Thus, \hat{x} will clearly be negative when 4\kappa\theta + Z^2 < 1. \blacksquare

Let Y(t) be such that d\log(Y(t)) = -\frac{1}{2}x(t)dt +\sqrt{x(t)}dW_Y(t) and dx(t) = \kappa(\theta - x(t))dt + \epsilon\sqrt(x(t))dW(t) (a CIR process).  Now define  A(T) = \frac{1}{T} \int_0^TY(s)ds.  We will (with our estimated parameters) use the antithetic method and (10) to estimate E(\max(3, A(T)) by simulation.

We wish to simulate 2 processes (Y_1(t)), (Y_2(t)) that are identically distributed and negatively correlated. To do this, we first simulate 2 processes (x_1(t)), (x_2(t)) that are identically distributed and negatively correlated using (10), the Milstein equation.

Next then obtain the processes d\log(Y_1(t)), d\log(Y_2(t)) and from there, (Y_1(t)), (Y_2(t)), trivially.

Recall that we wish to estimate E(\max(3, A(T)). We calculate \bar{A}(t)=\frac{A_1(t) +A_2(t)}{2}. Taking \max(3, A(T)) we obtain a sample value m_1. To estimate E(\max(3, A(T)), we repeat this simulation n times, and E(\max(3, A(T)) \approx \frac{1}{n}\sum_{i}m_i.

R Code

After running this experiment with our parameters \sigma = 0.98, \beta = -2.42, \gamma = .2, we obtain E(\max(3, A(T)) \approx 3.  Looking more closely at our data, we see that (indeed)  \bar{A_i}(t) < 3 \forall i.

Let us inspect a realization of \hat{x}_2(t):

This is what we expect (mean reversion at 0.2).  Let us now see the Heston process that was adapted from this realization.

 

Recall that (or rather, infer from d\log(Y(t))) dY(t) = Y(t)\sqrt{x(t)}dW_Y(t).  x(t) is a strictly positive process that hovers around its mean.  For the sake of illustration, we can (roughly) think of as dY(t) = Y(t)d\tilde{W}(t).  Then obviously, Y(t) is a log-normal process, so this plot is what we would expect.

In summary, despite our uninteresting result for E(\max(3, A(T)), we have shown how to estimate the parameters defining a CIR process from a (simulated CIR) dataset, and how to use antithetics to estimate an expectation dealing with the Heston model, where the underlying volatility is modeled by a CIR model.

This entry was posted in Uncategorized. Bookmark the permalink.

Leave a comment