added stochastic volatility example: source, optimized source, sim, and manual

Author: Bob Carpenter

misc/moving-avg/stochastic-volatility-optimized.stan

@@ -0,0 +1,25 @@
+
+data {
+  int<lower=0> T;   // # time points (equally spaced)
+  vector[T] y;      // mean corrected return at time t
+}
+parameters {
+  real mu;                     // mean log volatility
+  real<lower=-1,upper=1> phi;  // persistence of volatility
+  real<lower=0> sigma;         // white noise shock scale
+  vector[T] h_std;             // std log volatility time t
+}
+transformed parameters {
+  vector[T] h;                 // log volatility at time t
+  h <- h_std * sigma;
+  h[1] <- h[1] / sqrt(1 - phi * phi);
+  h <- h + mu;
+  for (t in 2:T)
+    h[t] <- h[t] + phi * (h[t-1] - mu);
+}
+model {
+  sigma ~ cauchy(0,5);
+  mu ~ cauchy(0,10);  
+  h_std ~ normal(0,1);
+  y ~ normal(0, exp(h / 2));
+}

misc/moving-avg/stochastic-volatility-sim.R

@@ -0,0 +1,20 @@
+
+phi <- 0.95;
+sigma <- 0.25;
+beta <- 0.6;
+mu <- 2 * log(beta);
+
+T <- 500;
+
+h <- rep(NA,T);
+h[1] <- rnorm(1, mu, sigma / sqrt(1 - phi * phi));
+for (t in 2:T)
+  h[t] <- rnorm(1, mu + phi * (h[t-1] - mu), sigma);
+y <- rep(NA,T);
+for (t in 1:T)
+  y[t] <- rnorm(1, 0, exp(h[t] / 2));
+
+library("rstan");
+start_time <- proc.time();
+fit <- stan(file="stochastic-volatility.stan", data=list(T=T,y=y), iter=10000, chains=4, init=0, fit=fit);
+elapsed_time <- proc.time() - start_time;

misc/moving-avg/stochastic-volatility.stan

@@ -0,0 +1,20 @@
+data {
+  int<lower=0> T;   // # time points (equally spaced)
+  vector[T] y;      // mean corrected return at time t
+}
+parameters {
+  real mu;                     // mean log volatility
+  real<lower=-1,upper=1> phi;  // persistence of volatility
+  real<lower=0> sigma;         // white noise shock scale
+  vector[T] h;                 // log volatility at time t
+}
+model {
+  phi ~ uniform(-1,1);
+  sigma ~ cauchy(0,5);
+  mu ~ cauchy(0,10);  
+  h[1] ~ normal(mu, sigma / sqrt(1 - phi * phi));
+  for (t in 2:T)
+    h[t] ~ normal(mu + phi * (h[t - 1] -  mu), sigma);
+  for (t in 1:T)
+    y[t] ~ normal(0, exp(h[t] / 2));
+}