Autoregressive Model
Simulation
\[\begin{align} Y_t &= x_t + v_t, \quad v_t \sim \mathcal{N}(0, V), \\ \alpha_t &= \phi (\alpha_{t-1}-\mu) + eta_t, \quad w_t \sim \mathcal{N}(0, \sigma_\eta), \\ \alpha_0 &\sim \mathcal{N}(m_0, C_0). \end{align}\]
The code required to simulate from this model is given below:
import com.github.jonnylaw.dlm._
import breeze.linalg.{DenseMatrix, DenseVector, diag}
val mod = Dlm.autoregressive(phi = 0.9)
val p = DlmParameters(
DenseMatrix(4.0),
DenseMatrix(2.0),
DenseVector(0.0),
DenseMatrix(1.0))
val data = Dlm.simulateRegular(mod, p, 1.0).
steps.
take(1000).
toVectorThe figure below shows a plot of 100 simulations:
Filtering
To perform Kalman Filtering on a Vector of Data, we simply discard the state from the simulated data and pass it into the KalmanFilter.filterDlm function:
import cats.implicits._
val filtered = KalmanFilter.filterDlm(mod, data.map(_._1), p)The filterDlm function is implemented for any collection which has a Traverse instance, cats provides an instance for Vector[A] (via the import cats.implicits._).
Parameter Inference
The figure below shows the prior distributions for the parameters and the diagnostics from 100,000 iterations of the MCMC algorithm. A Beta prior and proposal is used in a Metropolis-Hastings step to learn the autoregressive parameter \(\phi\).
import breeze.stats.distributions._
val priorV = InverseGamma(5.0, 20.0)
val priorW = InverseGamma(6.0, 10.0)
val priorPhi = new Beta(20, 2)First we define the prior distributions for the parameters, the prior for V and W are plotted in the Figure below
Next we can use the prior distribution to construct at DlmParameters object which can be drawn from to initialise the Markov chain.
val prior = for {
v <- priorV
w <- priorW
} yield DlmParameters(DenseMatrix(v), DenseMatrix(w), p.m0, p.c0)Next we create a single step of the Markov Chain which samples the value of \(\phi\) using a Metropolis Hastings step. The function updateModel is used to change the value of the system evolution matrix \(G\) to include the new value of the autoregressive parameter \(\phi\).
val step = (s: (Double, GibbsSampling.State)) => for {
newS <- GibbsSampling.dinvGammaStep(GibbsSampling.updateModel(mod, s._1),
priorV, priorW, data.map(_._1))(s._2)
phi <- GibbsSampling.samplePhi(priorPhi, 1000, 0.5, newS)(s._1)
} yield (phi, newS)Next the inital parameter distribution can be defined as:
val init = for {
p <- prior
phi <- priorPhi
state <- Smoothing.ffbsDlm(mod, data.map(_._1), p)
} yield (phi, GibbsSampling.State(p, state))This allows us to intialise the MCMC with a draw from the prior distribution. Finally we can construct the Markov Chain:
val iters = MarkovChain(init.draw)(step).steps.take(100000)The figure below shows the posterior distributions of the static parameters in the autoregressive DLM. 10,000 iterations are discarded as burnin and 100,000 samples are taken from the Markov Chain.
Continuous Time Autoregressive Process
The Ornstein-Uhlenbeck process is a continuous time autoregressive process. Define a partially observed Markov process with an OU latent-state as:
\[\begin{align*} Y_t &= x_t + v_t, \quad v_t \sim \mathcal{N}(0, V), \\ \textrm{d}\alpha_t &= \phi(\alpha_t - \mu)\textrm{d}t + \sigma \textrm{d}W_t. \end{align*}\]The Figure below shows a simulation from the OU process DLM with observation variance \(V = 0.5\), \(\phi = 0.2\), \(\mu = 1.0\) and \(\sigma = 0.3\). To simulate the OU process DLM, define parameters of the OU process and use the transition kernel defined in the StochasticVolatility class. The observation distribution is Gaussian with variance \(V = 0.5\). Since this process can be simulated at arbitrary time points, a vector of random times is simulated by defining the difference between the observations to be deltas. Then the state, \(\alpha_0\) is initialised at the stationary solution of the OU process and a scan is used to simulate the data.
scala> import com.github.jonnylaw.dlm._
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:22: warning: Unused import
import com.github.jonnylaw.dlm._
^
import com.github.jonnylaw.dlm._
scala> val p = SvParameters(0.2, 1.0, 0.3)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
p: com.github.jonnylaw.dlm.SvParameters = SvParameters(0.2,1.0,0.3)
scala> def stepDlm(t: Double, dt: Double, x: Double) =
| for {
| x1 <- StochasticVolatility.stepOu(p, x, dt)
| y <- Gaussian(x1, math.sqrt(0.5))
| } yield (t + dt, y, x1)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
stepDlm: (t: Double, dt: Double, x: Double)breeze.stats.distributions.Rand[(Double, Double, Double)]
scala> val deltas = Vector.fill(5000)(scala.util.Random.nextDouble())
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:20: warning: Unused import
import com.github.jonnylaw.dlm._
^
deltas: scala.collection.immutable.Vector[Double] = Vector(0.9112652524891296, 0.7170770448866899, 0.8861592330571312, 0.27463750868967585, 0.2616763853561149, 0.40735279441147565, 0.4115805846130032, 0.1714318530384753, 0.6422959068219808, 0.12274509299680303, 0.12991670749722417, 0.9703919078187349, 0.43335889663980975, 0.4898222138198741, 0.3131044366055894, 0.18608258678301481, 0.7254068467362147, 0.40424292348186364, 0.1463444483450903, 0.8277156979175803, 0.7237703499815168, 0.1651084945599559, 0.46367745798055904, 0.8971764949998614, 0.6770627693085174, 0.09771893596462333, 0.08095055010925845, 0.9318664758150106, 0.2403594414452308, 0.6137387001307241, 0.8248371072146025, 0.5301365992975726, 0.4425267131293412, 0.330628432879301, 0.5059874880342983, 0.9855625510887926, 0.7172039...
scala> val init = Gaussian(p.mu, math.sqrt(p.sigmaEta * p.sigmaEta / 2 * p.phi))
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:20: warning: Unused import
import com.github.jonnylaw.dlm._
^
init: breeze.stats.distributions.Gaussian = Gaussian(1.0, 0.09486832980505137)
scala> val sims = deltas.scanLeft((0.0, 0.0, init.draw)) {
| case ((t, y, xt), dt) =>
| stepDlm(t + dt, dt, xt).draw
| }
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:20: warning: Unused import
import com.github.jonnylaw.dlm._
^
sims: scala.collection.immutable.Vector[(Double, Double, Double)] = Vector((0.0,0.0,1.2674664723598252), (1.8225305049782592,0.7097903517586646,0.7363591030582783), (3.256684594751639,0.8982230407590618,0.9438613298682959), (5.029003060865901,-0.2956180673596438,0.9278377972558407), (5.578278078245252,0.9267980411369601,1.0300856169569532), (6.101630848957481,0.8765858414149623,0.7736046688534646), (6.916336437780433,0.06618320263200961,0.9607844474734939), (7.739497607006439,1.1523750169700686,0.8342216698015905), (8.08236131308339,1.247341169794061,0.9557244727919915), (9.36695312672735,1.4070041544885619,0.9151253216412034), (9.612443312720956,0.4407009735294981,0.9738279180733308), (9.872276727715406,0.42433406504140525,1.0164008947436398), (11.813060543352876,1.6205818934851994,1.0...
The figure below shows a Kalman Filter for the OU process, since the transition kernel of the OU process is known and given by a Gaussian distribution the Kalman Filter can be performed:
scala> val ys = sims.map { case (t, y, a) => (t, y.some) }
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:20: warning: Unused import
import com.github.jonnylaw.dlm._
^
ys: scala.collection.immutable.Vector[(Double, Option[Double])] = Vector((0.0,Some(0.0)), (1.8225305049782592,Some(0.7097903517586646)), (3.256684594751639,Some(0.8982230407590618)), (5.029003060865901,Some(-0.2956180673596438)), (5.578278078245252,Some(0.9267980411369601)), (6.101630848957481,Some(0.8765858414149623)), (6.916336437780433,Some(0.06618320263200961)), (7.739497607006439,Some(1.1523750169700686)), (8.08236131308339,Some(1.247341169794061)), (9.36695312672735,Some(1.4070041544885619)), (9.612443312720956,Some(0.4407009735294981)), (9.872276727715406,Some(0.42433406504140525)), (11.813060543352876,Some(1.6205818934851994)), (12.679778336632495,Some(0.7825244871108019)), (13.659422764272243,Some(-0.1743629107022191)), (14.285631637483423,Some(-0.7103486761484656)), (14.657796...
scala> val filtered = FilterOu.filterUnivariate(ys, Vector.fill(ys.size)(0.5), p)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
filtered: scala.collection.immutable.Vector[com.github.jonnylaw.dlm.FilterAr.FilterState] = Vector(FilterState(0.0,1.0,0.09,1.0,0.09), FilterState(0.0,0.847457627118644,0.07627118644067797,1.0,0.09), FilterState(1.8225305049782592,0.8508247205567557,0.11730129618998475,0.8940532815700306,0.1532554134913103), FilterState(3.256684594751639,0.8905461395477265,0.12367321421817269,0.8880232602738828,0.16431625237788855), FilterState(5.029003060865901,0.6070459365501255,0.12970084279855587,0.9232126892800705,0.1751298109598426), FilterState(5.578278078245252,0.7117853845806112,0.11449406893664624,0.647927290739849,0.14849845321553606), FilterState(6.101630848957481,0.769436654126097,0.10652637015317834,0.7404278130731669,0.13536659393750072), FilterState(6.916336437780433,0.6431569860212298,0...
Parameter inference is performed using the Kalman filter to calculate the marginal likelihood of the OU process which is then used in a Metropolis-Hastings algorithm to determine the static parameters of the latent state, \(\phi, \mu\) and \(\sigma\). The value of the observation variance, \(V\), is determined using a Gibbs step:
scala> import StochasticVolatility._
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:25: warning: Unused import
import StochasticVolatility._
^
import StochasticVolatility._
scala> val p = SvParameters(0.2, 1.0, 0.3)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
p: com.github.jonnylaw.dlm.SvParameters = SvParameters(0.2,1.0,0.3)
scala> val priorPhi = new Beta(2.0, 5.0)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
priorPhi: breeze.stats.distributions.Beta = breeze.stats.distributions.Beta@163831a8
scala> val priorMu = Gaussian(1.0, 1.0)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
priorMu: breeze.stats.distributions.Gaussian = Gaussian(1.0, 1.0)
scala> val priorSigma = InverseGamma(5.0, 1.0)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
priorSigma: com.github.jonnylaw.dlm.InverseGamma = InverseGamma(5.0,1.0)
scala> val priorV = InverseGamma(2.0, 2.0)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
priorV: com.github.jonnylaw.dlm.InverseGamma = InverseGamma(2.0,2.0)
scala> val f = (dt: Double) => DenseMatrix(1.0)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
f: Double => breeze.linalg.DenseMatrix[Double] = $$Lambda$6738/1042733082@23426493
scala> val step = (s: (StochasticVolatilityKnots.OuSvState, DenseMatrix[Double])) =>
| for {
| theta <- FilterOu.ffbs(s._1.params, ys, Vector.fill(ys.size)(s._2(0, 0)))
| st = theta.map(x => (x.time, x.sample))
| (phi, acceptedPhi) <- samplePhiOu(priorPhi, s._1.params, st, 0.05, 0.25)(
| s._1.params.phi)
| (mu, acceptedMu) <- sampleMuOu(priorMu, 0.2, s._1.params, st)(
| s._1.params.mu)
| (sigma, acceptedSigma) <- sampleSigmaMetropOu(priorSigma,
| 0.1,
| s._1.params,
| st)(s._1.params.sigmaEta)
| v <- GibbsSampling.sampleObservationMatrix(
| priorV,
| f,
| ys.map(x => DenseVector(x._2)),
| st.map { case (t, x) => (t, DenseVector(x)) })
| accepted = DenseVector(Array(acceptedPhi, acceptedMu, acceptedSigma))
| } yield
| (StochasticVolatilityKnots.OuSvState(SvParameters(phi, mu, sigma),
| theta,
| s._1.accepted + accepted),
| v)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
step: ((com.github.jonnylaw.dlm.StochasticVolatilityKnots.OuSvState, breeze.linalg.DenseMatrix[Double])) => breeze.stats.distributions.Rand[(com.github.jonnylaw.dlm.StochasticVolatilityKnots.OuSvState, breeze.linalg.DenseMatrix[Double])] = $$Lambda$6781/1487300472@77c16fae
scala> val initState = FilterOu.ffbs(p, ys, Vector.fill(ys.size)(priorV.draw))
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:24: warning: Unused import
import StochasticVolatility._
^
initState: breeze.stats.distributions.Rand[Vector[com.github.jonnylaw.dlm.FilterAr.SampleState]] = breeze.stats.distributions.RandBasis$$anon$9@78819b1
scala> val init = (StochasticVolatilityKnots.OuSvState(p,
| initState.draw,
| DenseVector.zeros[Int](3)),
| DenseMatrix(priorV.draw))
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:17: warning: Unused import
import breeze.stats.distributions._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
init: (com.github.jonnylaw.dlm.StochasticVolatilityKnots.OuSvState, breeze.linalg.DenseMatrix[Double]) = (OuSvState(SvParameters(0.2,1.0,0.3),Vector(SampleState(0.0,1.282651836274979,1.0,0.09,1.0,0.09), SampleState(0.0,1.282651836274979,0.9569170058287079,0.08612253052458371,1.0,0.09), SampleState(1.8225305049782592,1.4114428394695446,0.9372100455695462,0.13805552091956633,0.9700771545219369,0.1580075566096492), SampleState(3.256684594751639,1.1509840174489756,0.9484432680810332,0.16176010609461142,0.9528674294365996,0.1760103962350349), SampleState(5.029003060865901,0.7522288442732704,0.7414913596161166,0.15964900984398722,0.9638303959567323,0.19387510602841157), SampleState(5.578278078245252,0.9126935798009213,0.8014307919717945,0.13654724647161118,0.7683855548201457,0.172539390928050...
scala> val iters = MarkovChain(init)(step)
<console>:10: warning: Unused import
import com.github.jonnylaw.dlm._
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:12: warning: Unused import
import breeze.linalg.{DenseMatrix, DenseVector, diag}
^
<console>:14: warning: Unused import
import cats.implicits._
^
<console>:23: warning: Unused import
import StochasticVolatility._
^
iters: breeze.stats.distributions.Process[(com.github.jonnylaw.dlm.StochasticVolatilityKnots.OuSvState, breeze.linalg.DenseMatrix[Double])] = breeze.stats.distributions.MarkovChain$$anon$1@734b0025