In [1]:
import Resolvers._
interp.resolvers() = interp.resolvers() :+ Resolver.Http(
"cibotech",
"https://dl.bintray.com/cibotech/public/",
MavenPattern,
true
)
Out[1]:
In [2]:
import $ivy.`com.cibo::evilplot:0.4.1`
import $ivy.`com.cibo::evilplot-repl:0.4.1`
import $ivy.`org.ruivieira::scala-ssm:0.1`
Out[2]:
In [3]:
import breeze.linalg.{DenseVector, DenseMatrix}
import org.ruivieira.ssm.univariate.{
StateGenerator,
UnivariateGenerator,
UnivariateStructure
}
Out[3]:
First, we start by creating a simple univariate Gaussian random walk. This will correspond to a dynamic linear model in the form
\begin{align} y_t &\sim \mathcal{N}\left(\theta_t,V\right) \\ \theta_t &\sim \mathcal{N}\left(\theta_{t-1},W\right) \end{align}We start by definind the variance of the latent states as $W=1.5$.
In [4]:
val structure = UnivariateStructure.createLocallyConstant(W = 1.5)
Out[4]:
And generate a chain of $n=1000$ states with an initial state of $\theta_0 = 0$.
In [5]:
val states = StateGenerator.states(nobs = 1000,
structure = structure,
state0 = DenseVector[Double](0.0))
Out[5]:
We can now generate the observations from the states, using an observation variance $V=4.0$.
In [10]:
val observations = UnivariateGenerator.gaussian(states = states,
structure = structure,
V = 4.0)
Out[10]:
In [11]:
import com.cibo.evilplot._
import com.cibo.evilplot.plot._
import com.cibo.evilplot.plot.aesthetics.DefaultTheme._
import com.cibo.evilplot.numeric.Point
val states_plot = ScatterPlot(Seq.tabulate(100) { i => Point(i.toDouble, states(i)(0)) })
val obs_plot = LinePlot(Seq.tabulate(100) { i => Point(i.toDouble, observations(i)) })
val plot = Overlay(states_plot, obs_plot).render()
publish.png(plot.asBufferedImage)
Out[11]:
Let's assume we are following the end of day stock price of company Foo Ltd. We will simulate a stock price history for 365 days (each timepoint is a day), for a relatively stable stock price, with no trend or seasonality but with natural fluctuation. We also assume that the stock's initial price at day 0 is around \$100. As such we will set an initial state of $\theta_0 = 100$ and low underlying variance ($W=0.01$) with some noise in the data ($V=1.0$).
In [7]:
val stock_structure = UnivariateStructure.createLocallyConstant(W = 0.01)
val states = StateGenerator.states(nobs = 365,
structure = structure,
state0 = DenseVector[Double](100.0))
val observations = UnivariateGenerator.gaussian(states = states,
structure = structure,
V = 1.0)
Out[7]:
In [8]:
Scatter((1 until 365), observations.toSeq).plot()
Out[8]:
Now (ignoring the complexities of the stock market), let's suppose that on day $t=100$, this company announces a revolutionary breakthrough. We want to incorporate in our simulated data a jump of twice is stock price, regardless of the value at $t=100$.
We can do this by changing the data at the state level at any point. First we create a chain for a "normal" random walk and then we append another one with a starting value of $\theta_{100} = 2\theta_{100}$.
In [9]:
val states_pre = StateGenerator.states(nobs = 100, structure = structure, state0 = DenseVector[Double](100.0))
val states_post = StateGenerator.states(nobs = 265, structure = structure, state0 = states_pre.last * 2.0)
val states = states_pre ++ states_post
Out[9]:
It is important to note that due to the Markovian nature of the SSM, changes at the state level will propagate to future states. This means that the jump in stock price will be propagated to future values.
In [10]:
val observations = UnivariateGenerator.gaussian(states = states,
structure = structure,
V = 1.0)
Scatter((1 until 365), observations.toSeq).plot()
Out[10]:
For a locally linear model, we assume the state and observations matrices to be, respectively
$$ \mathsf{F} = \begin{bmatrix} 1 & 0 \end{bmatrix},\qquad \mathsf{G} = \begin{bmatrix} 1 & 1 \\ 0 & 1\end{bmatrix}.$$The latent states, will then correspond to $\theta_t = \left(\mu, \tau\right)$, that is two components, representing the mean and the trend, respectively. The model will then take the form
\begin{align} y_t \sim \mathcal{N}\left(\mathsf{F}\theta_t,V\right) \\ \theta_t \sim \mathcal{N}\left(\mathsf{G}\theta_{t-1},\mathsf{W}\right) \end{align}The state covariance will now be a matrix
$$ \mathsf{W} = \begin{bmatrix} W_{\tau} & 0 \\ 0 & W_{\mu} \end{bmatrix}$$representing the variance of the underlying mean and trend respectively.
Let's now simulate the stock price of the Foo company, but assuming there's a trend to the values, rather than a jump. We will create a mean varying a bit ($W_{\mu}=0.5$) but a rather smooth trend ($W_{\tau}=0.05$) and with some noise in the observations ($V=2.0$).
In [27]:
val W = DenseMatrix.eye[Double](2)
W(0,0) = 0.05
W(1,1) = 0.5
val stock_structure = UnivariateStructure.createLocallyLinear(W = W)
val states = StateGenerator.states(nobs = 365,
structure = stock_structure,
state0 = DenseVector[Double](100.0, -1.0))
val observations = UnivariateGenerator.gaussian(states = states,
structure = stock_structure,
V = 2.0)
Out[27]:
In [28]:
Scatter((1 until 365), observations.toSeq).plot()
Out[28]:
One of the advantages of this formulation, is that we can decompose easily the states into the mean and the trend. For instance:
In [31]:
val x = 1 until 365
val plot = Seq(
Scatter(
x, states.map(_(0)), name = "trend"
),
Scatter(
x, states.map(_(1)), name = "mean"
)
)
plot.plot(title = "Locally linear")
Out[31]:
In [ ]: