User:Carl Boettiger/Notebook/Stochastic Population Dynamics/2010/05/06

{| width="800"
 * style="background-color: #EEE"|[[Image:owwnotebook_icon.png|128px]] Stochastic Population Dynamics
 * style="background-color: #F2F2F2" align="center"|  |Main project page
 * style="background-color: #F2F2F2" align="center"|  |Main project page


 * colspan="2"|
 * colspan="2"|

Goals

 * Model choice exercise with sdes
 * general likelihood calculation from simulation

Test Case
$$ dX_t = \alpha_t (\theta - X_t) dt + \sigma dt $$

$$ d \alpha_t = -\beta $$

Approach

 * The linear approximation to the warning signals dynamics can be captured by the OU process, look for changing $$\alpha$$
 * SDE models will also provide some of the coarser approximations for the structured population dynamics. Formulations of these still to do.

Conditional probability
Solution to the time dependent OU process, see Gardiner 4th ed. pg 111, eq. (4.5.109).

Solve using $$ f(X_t) = X_t \exp\left( - \int_0^t \alpha(s) ds \right) $$ and apply Ito formula, giving


 * $$ X_t = X_0 \exp\left( - \int_0^t \alpha(s) ds \right) +\theta\left( 1- \exp\left( - \int_0^t \alpha(s) ds \right) \right) + \int_0^t \sigma \exp\left( - \int_0^t \alpha(s) ds \right) dW_t $$

And the moments are:



\begin{align} E(X_T) &= X_0 \exp\left( - \int_0^T \alpha(s) ds \right) +\theta\left( 1- \exp\left( - \int_0^T \alpha(s) ds \right) \right) \\ \operatorname{Var}(X_T) &= \sigma^2 \int_0^T \exp\left( -2 \int_0^t \alpha(s) ds \right) d t \end{align} $$


 * Using the exact solution I can write down the conditional density function. We can evaluate both integrals for this linear increase model:


 * $$ \int_0^t (\beta t +\alpha_0 ) dt = \beta t^2/2 + \alpha_0 t $$


 * $$ \int_0^t e^{ - \beta s^2 - 2\alpha_0 s} ds = \frac{\sqrt{\pi}e^{\alpha_0^2/\beta} }{2\sqrt{\beta}} \left( \operatorname{erf} \left( \frac{ \alpha_0+\beta t }{ \sqrt{\beta} } \right) - \operatorname{erf} \left( \frac{ \alpha_0 }{ \sqrt{\beta} } \right) \right) $$


 * Next, implement the conditional probability density in R:

Note on Numerics

 * Unfortunately the analytic solution does not have nice numerical properties for small β, which is the interesting limiting case (no/very slow change in stability). The error functions both approach unity, making the difference tiny, while the exponential term diverges.  The result is that variances for too small beta suddenly become zero, instead of the appropriate limit.  The numerical integration is more robust to this.  Analytic implementation is modified to just provide the limit case when erf functions both evaluate at exactly 1.


 * This linear test case model is quite restrictive. If the eigenvalue is decreasing, doing so in a linear way need not out perform a model in which it is not decreasing at all.  It would be interesting to test if this is more sensitive than change point estimation.  Iacus (2008) implements a least-squares based change point estimation in the R sde package.

Misc Notes

 * Upgraded to Ubuntu 10.04. Went fine on the experimental partition with a reformat to extension 4 and a fresh install.  Fast boot times, spiffy system.  Full upgrade on main (ext3) partition failed miserably though, can't even boot up successfully following the upgrade.


 * Adding additional packages when fine, mostly copied configuration files over. A couple extra steps to get dual monitors working again, as the resolution wasn't properly detected.  Solution here, essentially:


 * vim-full is deprecated, vim-latexsuite needs additional installation, read: /usr/share/doc/vim-latexsuite/README.Debian after installing. Also vim's "+y copy to clipboard doesn't work with pure vim, so installing xclip vim-common and vim-gnome to fixed this.


 * }