In [9]:

```
import numpy as np
import pymc as pm
from pymc.distributions.timeseries import GaussianRandomWalk
from scipy.sparse import csc_matrix
from scipy import optimize
%pylab inline
```

Asset prices have time-varying volatility (variance of day over day `returns`

). In some periods, returns are highly variable, while in others very stable. Stochastic volatility models model this with a latent volatility variable, modeled as a stochastic process. The following model is similar to the one described in the No-U-Turn Sampler paper, Hoffman (2011) p21.

$$ \sigma \sim Exponential(50) $$

$$ \nu \sim Exponential(.1) $$

$$ s_i \sim Normal(s_{i-1}, \sigma^{-2}) $$

$$ log(\frac{y_i}{y_{i-1}}) \sim t(\nu, 0, exp(-2 s_i)) $$

Here, $y$ is the daily return series and $s$ is the latent log volatility process.

First we load some daily returns of the S&P 500.

In [12]:

```
n = 400
returns = np.genfromtxt("data/SP500.csv")[-n:]
returns[:5]
```

Out[12]:

In [24]:

```
plt.plot(returns)
```

Out[24]:

Specifying the model in pymc mirrors its statistical specification.

However, it is easier to sample the scale of the log volatility process innovations, $\sigma$, on a log scale, so we create it using `TransformedVar`

and use `logtransform`

. `TransformedVar`

creates one variable in the transformed space and one in the normal space. The one in the transformed space (here $\text{log}(\sigma) $) is the one over which sampling will occur, and the one in the normal space is the one to use throughout the rest of the model.

It takes a variable name, a distribution and a transformation to use.

In [20]:

```
model = pm.Model()
with model:
sigma, log_sigma = model.TransformedVar('sigma', pm.Exponential.dist(1./.02, testval=.1),
pm.logtransform)
nu = pm.Exponential('nu', 1./10)
s = GaussianRandomWalk('s', sigma**-2, shape=n)
r = pm.T('r', nu, lam=pm.exp(-2*s), observed=returns)
```

For this model, the full maximum a posteriori (MAP) point is degenerate and has infinite density. However, if we fix `log_sigma`

and `nu`

it is no longer degenerate, so we find the MAP with respect to the volatility process, 's', keeping `log_sigma`

and `nu`

constant at their default values.

We use L-BFGS because it is more efficient for high dimensional functions (`s`

has n elements).

In [21]:

```
with model:
start = pm.find_MAP(vars=[s], fmin=optimize.fmin_l_bfgs_b)
```

In [33]:

```
with model:
step = pm.NUTS(scaling=start)
start2 = pm.sample(500, step, progressbar=False)[-1]
# Start next run at the last sampled position.
step = pm.NUTS(scaling=start2)
trace = pm.sample(2000, step, start=start2, progressbar=False)
```

In [34]:

```
figsize(12,6)
pm.traceplot(trace, model.vars[:-1]);
```

In [35]:

```
figsize(12,6)
title(str(s))
plot(trace[s][::10].T,'b', alpha=.03);
xlabel('time')
ylabel('log volatility')
```

Out[35]:

In [36]:

```
plot(returns)
plot(np.exp(trace[s][::10].T), 'r', alpha=.03);
sd = np.exp(trace[s].T)
plot(-np.exp(trace[s][::10].T), 'r', alpha=.03);
xlabel('time')
ylabel('returns')
```

Out[36]:

- Hoffman & Gelman. (2011). The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.