# DynamicDiscreteModels.jl

Back-end support for doing statistical inference with any model that can described as a "partially observed" (discrete) Markov chain, such as a Hidden Markov Model.

Model construction and data simulation:

```
julia> #use a constructor defined in a front-end package:
julia> model=myawesomemodel()
julia> coef!(model,trueparameter)
julia> #or simply:
julia> model=myawesomemodel(trueparameter)
julia> #generate data from the model:
julia> data=rand(model,1000); #1 time series with 1000 observations
julia> data=rand(model,100,10); #10 time series with 100 observations each
julia> data=rand(model,[100,120,115]); #3 time series with 100, 120 and 115 obs.
```

Parameter estimation:

```
julia> #plain maximum likelihood:
julia> thetahat=mle(model,data)
julia> #maximum likelihood via EM algorithm:
julia> thetahat=em(model,data)
```

Hidden state inference (Viterbi filtering):

```
julia> #pretend thetahat is the "true" parameter:
julia> coef!(model,thetahat)
julia> #most likely path of latent states:
julia> viterbi(model,data)
```

## Scope of the package and role of front-end packages

The scope of the package includes:

- filter/smoother algorithms (aka. forward/backward)
- parameter estimation: maximum likelihood and EM algorithm
- hidden state inference: Viterbi filtering
- support for deep modelling of the transition matrices
- support for sparse transition matrices
- support for Jacobians in all of the above

In practice DynamicDiscreteModels.jl implements:

- some core methods:
`rand()`

,`loglikelihood()`

,`estep()`

and`viterbi()`

. - convenience methods:
`em()`

wraps`estep()`

with generic M-step numerical optimizations. DynamicDiscreteModels.jl also inherits`mle()`

from`StatisticalModels`

, which wraps`loglikelihood()`

with generic numerical optimization.

The job of a front-end package MyModel.jl is simply to specify `coef!(mymodel,myparameter)`

which maps `myparameter`

to DynamicDiscreteModels.jl's canonical representation. All of the above methods will then be available. Optionally, a front-end package may also extend `em()`

and/or `mle()`

with specialized optimization wrapping DynamicDiscreteModels.jl's `estep()`

and `loglikelihood()`

.

See **Usage** below for a simple example and see HiddenMarkovModels.jl for an actual front-end package built on top of DynamicDiscreteModels.jl.

## Installation

`julia> Pkg.add("DynamicDiscreteModels")`

## Usage

For the purpose of this package, a "dynamic discrete model" is a discrete Markov chain z=(x,y), where x=1...dx is latent and y=1...dy is observed. The package defines an abstract type `DynamicDiscreteModel`

with two fields:

`m`

is a (dx,dy,dx,dy) array that holds the Markov probabilities m[x,y,x',y'] of moving from (x,y) today to (x',y') tomorrow.`mu`

is a (dx,dy) array that holds the joint initial distribution of the chain.

Any front-end implementation of a `DynamicDiscreteModel`

will specify a mapping from a statistical parameter θ to the transition matrix `m`

in a `coef!(model,parameter)`

function. Being a back-end package, `DynamicDiscreteModels.jl`

is agnostic as to the specific mapping, but examples/toymodel.jl provides a simple example. In this example (x,y) is a Hidden Markov model where x moves from today to tomorrow according to the Markov transition matrix A(θ) and y is drawn conditional on x according to the emission/transition matrix B(θ):

```
.
[ θ₁ 1-θ₁ ] [θ₂ (1-θ₂)/2 (1-θ₂)/2 ]
A(θ) = [ 1-θ₁ θ₁ ] B(θ) = [(1-θ₂)/2 (1-θ₂)/2 θ₂ ]
```

The statistical parameter θ is of dimension 2, x can take values 1, 2 and y can take values 1, 2, 3. Here is the code defining a toy model:

```
importall DynamicDiscreteModels
type ToyModel <: DynamicDiscreteModel
#DynamicDiscreteModel fields
m::Array{Float64,4} #the transition matrix given as m[x,y,x',y']
mu::Array{Float64,2} #initial distribution (dx,dy)
#DynamicDiscreteModel's technical fields
rho::Array{Float64,1}
phi::Array{Float64,1}
psi::Array{Float64,1}
end
dx,dy=2,3
toymodel()=ToyModel(
Array(Float64,dx,dy,dx,dy),
fill(1/6,2,3),
Array(Float64,1),
Array(Float64,dx),
Array(Float64,dx))
function coef!(model::ToyModel,theta::Tuple)
p1,p2=theta[1],theta[2]
a=[p1 1-p1;1-p1 p1]
p3=(1-p2)/2
b=[p2 p3 p3; p3 p3 p2]
#a convenience function is provided for calibrating models with hidden Markov structure
hmm2ddm!(model,a,b)
end
```

We can initiate a toy model, calibrate it to a true parameter value θ*=(.65,.5), and simulate 100 time series of length 100:

```
julia> thetastar=(.65, .5);
julia> model=toymodel();
julia> coef!(model,thetastar);
julia> data=simulate(model,100,100);
julia> data[1]
100-element Array{Int64,1}:
1
3
1
1
2
3
...
```

With the model and some data, we can find the most likely parameter value behind the data, according to the model, by computing the maximum likelihood estimator. Two methods are available: direct numerical optimization of the log-likelihood via `mle()`

and the EM algorithm via `em()`

. To help the optimizer, it is a good idea to first reparametrize the model on R rather than on [0,1]:

```
theta2eta(theta::Tuple)=[log(theta[1]/(1-theta[1])),log(theta[2]/(1-theta[2]))]
eta2theta(eta::Array)=(exp(eta[1])/(1+exp(eta[1])),exp(eta[2])/(1+exp(eta[2])))
coef!(model::ToyModel,eta::Array)=calibrate!(model,eta2theta(eta))
```

Notice how `coef!()`

will use multiple dispatch to choose between the `eta`

and `theta`

parametrizations. `mle()`

and `em()`

will dispatch on `parameter::Array`

by default so keep this signature for a parametrization which plays nicely with numerical optimization.

We then simply call:

```
julia> thetahat=eta2theta(mle(model,data))
(0.6762467624375702,0.5026000995341361)
julia> thetahat2=eta2theta(em(model,data))
(0.6763156818261421,0.5025930701695173)
```

Good news, up to some numerical precision error, `mle()`

and `em()`

agree on the maximum of the likelihood function. It also turns out that with this sample size, we can pinpoint the true parameter value (0.65,0.5) rather well. This is not the case with a smaller sample size:

```
julia> thetahat3=eta2theta(mle(model,data[1:10]))
(0.8420953672648955,0.4391327611992014)
```

Here is a plot of the likelihood surface using 10 and 100 time series, along with their maximums `thetahat`

and `thetahat3`

computed above (in blue) and the true parameter value (in green):

(See examples/toymodel.jl for the code used generate the picture.)