Skip to content

Instantly share code, notes, and snippets.

@ChadFulton

ChadFulton/_dStatespace.html

Last active August 29, 2015 14:03
Show Gist options
  • Save ChadFulton/92919e92be57b7855067 to your computer and use it in GitHub Desktop.
Save ChadFulton/92919e92be57b7855067 to your computer and use it in GitHub Desktop.
Download ZIP
Literate Programming Output: _statespace
Raw
_dStatespace.html
<!DOCTYPE html>
<html>
<head>
<meta http-equiv="content-type" content="text/html;charset=utf-8">
<title>_statespace.pyx</title>
<style>
/*--------------------- Layout and Typography ----------------------------*/
body {
font-family: 'Palatino Linotype', 'Book Antiqua', Palatino, FreeSerif, serif;
font-size: 16px;
line-height: 24px;
color: #252519;
margin: 0; padding: 0;
background: #f5f5ff;
}
a {
color: #261a3b;
}
a:visited {
color: #261a3b;
}
p {
margin: 0 0 15px 0;
}
h1, h2, h3, h4, h5, h6 {
margin: 40px 0 15px 0;
}
h2, h3, h4, h5, h6 {
margin-top: 0;
}
#container {
background: white;
}
#container, div.section {
position: relative;
}
#background {
position: absolute;
top: 0; left: 580px; right: 0; bottom: 0;
background: #f5f5ff;
border-left: 1px solid #e5e5ee;
z-index: 0;
}
#jump_to, #jump_page {
background: white;
-webkit-box-shadow: 0 0 25px #777; -moz-box-shadow: 0 0 25px #777;
-webkit-border-bottom-left-radius: 5px; -moz-border-radius-bottomleft: 5px;
font: 10px Arial;
text-transform: uppercase;
cursor: pointer;
text-align: right;
}
#jump_to, #jump_wrapper {
position: fixed;
right: 0; top: 0;
padding: 5px 10px;
}
#jump_wrapper {
padding: 0;
display: none;
}
#jump_to:hover #jump_wrapper {
display: block;
}
#jump_page {
padding: 5px 0 3px;
margin: 0 0 25px 25px;
}
#jump_page .source {
display: block;
padding: 5px 10px;
text-decoration: none;
border-top: 1px solid #eee;
}
#jump_page .source:hover {
background: #f5f5ff;
}
#jump_page .source:first-child {
}
div.docs {
float: left;
max-width: 500px;
min-width: 500px;
min-height: 5px;
padding: 10px 25px 1px 50px;
vertical-align: top;
text-align: left;
}
.docs pre {
margin: 15px 0 15px;
padding-left: 15px;
}
.docs p tt, .docs p code {
background: #f8f8ff;
border: 1px solid #dedede;
font-size: 12px;
padding: 0 0.2em;
}
.octowrap {
position: relative;
}
.octothorpe {
font: 12px Arial;
text-decoration: none;
color: #454545;
position: absolute;
top: 3px; left: -20px;
padding: 1px 2px;
opacity: 0;
-webkit-transition: opacity 0.2s linear;
}
div.docs:hover .octothorpe {
opacity: 1;
}
div.code {
margin-left: 580px;
padding: 14px 15px 16px 50px;
vertical-align: top;
}
.code pre, .docs p code {
font-size: 12px;
}
pre, tt, code {
line-height: 18px;
font-family: Monaco, Consolas, "Lucida Console", monospace;
margin: 0; padding: 0;
}
div.clearall {
clear: both;
}
/*---------------------- Syntax Highlighting -----------------------------*/
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
body .hll { background-color: #ffffcc }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #954121 } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #808080 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0040D0 } /* Generic.Traceback */
body .kc { color: #954121 } /* Keyword.Constant */
body .kd { color: #954121; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #954121; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #954121 } /* Keyword.Pseudo */
body .kr { color: #954121; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #219161 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #954121 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #954121; font-weight: bold } /* Name.Tag */
body .nv { color: #19469D } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #219161 } /* Literal.String.Backtick */
body .sc { color: #219161 } /* Literal.String.Char */
body .sd { color: #219161; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #219161 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #219161 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #954121 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #219161 } /* Literal.String.Single */
body .ss { color: #19469D } /* Literal.String.Symbol */
body .bp { color: #954121 } /* Name.Builtin.Pseudo */
body .vc { color: #19469D } /* Name.Variable.Class */
body .vg { color: #19469D } /* Name.Variable.Global */
body .vi { color: #19469D } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */
</style>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
tex2jax: {
inlineMath: [ ['$','$'] ],
displayMath: [ ['$$','$$'] ],
processEscapes: true
},
});
</script>
<script type="text/javascript"
src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</head>
<body>
<div id='container'>
<div id="background"></div>
<div class='section'>
<div class='docs'><h1>_statespace.pyx</h1></div>
</div>
<div class='clearall'>
<div class='section' id='section-0'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-0'>#</a>
</div>
<p>cython: boundscheck=False
cython: wraparound=False
cython: cdivision=False</p>
<p>State Space Models</p>
<p>Author: Chad Fulton<br />
License: Simplified-BSD</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-1'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-1'>#</a>
</div>
<h2>Constants</h2>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-2'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-2'>#</a>
</div>
<h3>Filters</h3>
<p>TODO note that only the conventional filter is implemented</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_CONVENTIONAL</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x01</span> <span class="c"># Durbin and Koopman (2012), Chapter 4</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_EXACT_INITIAL</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x02</span> <span class="c"># ibid., Chapter 5.6</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_AUGMENTED</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x04</span> <span class="c"># ibid., Chapter 5.7</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_SQUARE_ROOT</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x08</span> <span class="c"># ibid., Chapter 6.3</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_UNIVARIATE</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x10</span> <span class="c"># ibid., Chapter 6.4</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_COLLAPSED</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x20</span> <span class="c"># ibid., Chapter 6.5</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_EXTENDED</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x40</span> <span class="c"># ibid., Chapter 10.2</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FILTER_UNSCENTED</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x80</span> <span class="c"># ibid., Chapter 10.3</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-3'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-3'>#</a>
</div>
<h3>Inversion methods</h3>
<p>Methods by which the terms using the inverse of the forecast error
covariance matrix are solved.</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="kt">int</span> <span class="nf">INVERT_UNIVARIATE</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x01</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">SOLVE_LU</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x02</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">INVERT_LU</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x04</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">SOLVE_CHOLESKY</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x08</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">INVERT_CHOLESKY</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x10</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-4'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-4'>#</a>
</div>
<h3>Numerical Stability methods</h3>
<p>Methods to improve numerical stability</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="kt">int</span> <span class="nf">STABILITY_FORCE_SYMMETRY</span> <span class="o">=</span> <span class="mf">0</span><span class="n">x01</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-5'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-5'>#</a>
</div>
<p>Typical imports</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="k">import</span> <span class="nn">warnings</span>
<span class="k">cimport</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="k">cimport</span> <span class="nn">cython</span>
<span class="n">np</span><span class="o">.</span><span class="n">import_array</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-6'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-6'>#</a>
</div>
<h2>Math Functions</h2>
<p>Real and complex log and abs functions</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">from</span> <span class="nn">libc.math</span> <span class="k">cimport</span> <span class="n">log</span> <span class="k">as</span> <span class="n">dlog</span><span class="p">,</span> <span class="nb">abs</span> <span class="k">as</span> <span class="n">dabs</span>
<span class="k">from</span> <span class="nn">numpy</span> <span class="k">cimport</span> <span class="n">npy_cdouble</span>
<span class="k">cdef</span> <span class="kr">extern</span> <span class="k">from</span> <span class="s">&quot;numpy/npy_math.h&quot;</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">NPY_PI</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">npy_cabs</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">npy_cdouble</span> <span class="n">z</span><span class="p">)</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_cdouble</span> <span class="n">npy_clog</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">npy_cdouble</span> <span class="n">z</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kr">inline</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">zabs</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">complex128_t</span> <span class="n">z</span><span class="p">):</span>
<span class="k">return</span> <span class="n">npy_cabs</span><span class="p">((</span><span class="o">&lt;</span><span class="n">np</span><span class="o">.</span><span class="n">npy_cdouble</span> <span class="o">*&gt;</span> <span class="o">&amp;</span><span class="n">z</span><span class="p">)[</span><span class="mf">0</span><span class="p">])</span>
<span class="k">cdef</span> <span class="kr">inline</span> <span class="kt">np</span>.<span class="kt">complex128_t</span> <span class="nf">zlog</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">complex128_t</span> <span class="n">z</span><span class="p">):</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">npy_cdouble</span> <span class="nf">x</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">npy_clog</span><span class="p">((</span><span class="o">&lt;</span><span class="n">np</span><span class="o">.</span><span class="n">npy_cdouble</span><span class="o">*&gt;</span> <span class="o">&amp;</span><span class="n">z</span><span class="p">)[</span><span class="mf">0</span><span class="p">])</span>
<span class="k">return</span> <span class="p">(</span><span class="o">&lt;</span><span class="n">np</span><span class="o">.</span><span class="n">complex128_t</span> <span class="o">*&gt;</span> <span class="o">&amp;</span><span class="n">x</span><span class="p">)[</span><span class="mf">0</span><span class="p">]</span>
<span class="k">cdef</span> <span class="kr">extern</span> <span class="k">from</span> <span class="s">&quot;capsule.h&quot;</span><span class="p">:</span>
<span class="n">void</span> <span class="o">*</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="nb">object</span> <span class="n">ptr</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-7'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-7'>#</a>
</div>
<h2>BLAS / LAPACK functions</h2>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-8'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-8'>#</a>
</div>
<p><code>blas_lapack.pxd</code> contains typedef statements for BLAS and LAPACK functions</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">from</span> <span class="nn">blas_lapack</span> <span class="k">cimport</span> <span class="o">*</span>
<span class="k">try</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-9'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-9'>#</a>
</div>
<p>Scipy &gt;= 0.12.0 exposes Fortran BLAS functions directly</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">from</span> <span class="nn">scipy.linalg.blas</span> <span class="k">import</span> <span class="n">cgerc</span>
<span class="k">except</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-10'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-10'>#</a>
</div>
<p>Scipy &lt; 0.12.0 exposes Fortran BLAS functions in the <code>fblas</code> submodule</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">from</span> <span class="nn">scipy.linalg.blas</span> <span class="k">import</span> <span class="n">fblas</span> <span class="k">as</span> <span class="n">blas</span>
<span class="k">else</span><span class="p">:</span>
<span class="k">from</span> <span class="nn">scipy.linalg</span> <span class="k">import</span> <span class="n">blas</span>
<span class="k">try</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-11'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-11'>#</a>
</div>
<p>Scipy &gt;= 0.12.0 exposes Fortran LAPACK functions directly</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">from</span> <span class="nn">scipy.linalg.lapack</span> <span class="k">import</span> <span class="n">cgbsv</span>
<span class="k">except</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-12'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-12'>#</a>
</div>
<p>Scipy &lt; 0.12.0 exposes Fortran LAPACK functions in the <code>flapack</code> submodule</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">from</span> <span class="nn">scipy.linalg.lapack</span> <span class="k">import</span> <span class="n">flapack</span> <span class="k">as</span> <span class="n">lapack</span>
<span class="k">else</span><span class="p">:</span>
<span class="k">from</span> <span class="nn">scipy.linalg</span> <span class="k">import</span> <span class="n">lapack</span>
<span class="k">cdef</span> <span class="kt">dgemm_t</span> *<span class="nf">dgemm</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dgemm_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">dgemm</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dgemv_t</span> *<span class="nf">dgemv</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dgemv_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">dgemv</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dcopy_t</span> *<span class="nf">dcopy</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dcopy_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">dcopy</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">daxpy_t</span> *<span class="nf">daxpy</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">daxpy_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">daxpy</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dscal_t</span> *<span class="nf">dscal</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dscal_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">dscal</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dgetrf_t</span> *<span class="nf">dgetrf</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dgetrf_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dgetrf</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dgetri_t</span> *<span class="nf">dgetri</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dgetri_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dgetri</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dgetrs_t</span> *<span class="nf">dgetrs</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dgetrs_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dgetrs</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dpotrf_t</span> *<span class="nf">dpotrf</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dpotrf_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dpotrf</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dpotri_t</span> *<span class="nf">dpotri</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dpotri_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dpotri</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">dpotrs_t</span> *<span class="nf">dpotrs</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">dpotrs_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">lapack</span><span class="o">.</span><span class="n">dpotrs</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">ddot_t</span> *<span class="nf">ddot</span> <span class="o">=</span> <span class="o">&lt;</span><span class="n">ddot_t</span><span class="o">*&gt;</span><span class="n">Capsule_AsVoidPtr</span><span class="p">(</span><span class="n">blas</span><span class="o">.</span><span class="n">ddot</span><span class="o">.</span><span class="n">_cpointer</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">FORTRAN</span> <span class="o">=</span> <span class="mf">1</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-13'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-13'>#</a>
</div>
<p>Array shape validation</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="nf">validate_matrix_shape</span><span class="p">(</span><span class="nb">str</span> <span class="n">name</span><span class="p">,</span> <span class="n">Py_ssize_t</span> <span class="o">*</span><span class="n">shape</span><span class="p">,</span> <span class="nb">int</span> <span class="n">nrows</span><span class="p">,</span> <span class="nb">int</span> <span class="n">ncols</span><span class="p">,</span> <span class="n">nobs</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">==</span> <span class="n">nrows</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Invalid shape for </span><span class="si">%s</span><span class="s"> matrix: requires </span><span class="si">%d</span><span class="s"> rows,&#39;</span>
<span class="s">&#39; got </span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">name</span><span class="p">,</span> <span class="n">nrows</span><span class="p">,</span> <span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]))</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">==</span> <span class="n">ncols</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Invalid shape for </span><span class="si">%s</span><span class="s"> matrix: requires </span><span class="si">%d</span><span class="s"> columns,&#39;</span>
<span class="s">&#39;got </span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">name</span><span class="p">,</span> <span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">],</span> <span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]))</span>
<span class="k">if</span> <span class="n">nobs</span> <span class="ow">is</span> <span class="ow">not</span> <span class="bp">None</span> <span class="ow">and</span> <span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="ow">not</span> <span class="ow">in</span> <span class="p">[</span><span class="mf">1</span><span class="p">,</span> <span class="n">nobs</span><span class="p">]:</span>
<span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Invalid time-varying dimension for </span><span class="si">%s</span><span class="s"> matrix:&#39;</span>
<span class="s">&#39; requires 1 or </span><span class="si">%d</span><span class="s">, got </span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">name</span><span class="p">,</span> <span class="n">nobs</span><span class="p">,</span> <span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]))</span>
<span class="k">cdef</span> <span class="nf">validate_vector_shape</span><span class="p">(</span><span class="nb">str</span> <span class="n">name</span><span class="p">,</span> <span class="n">Py_ssize_t</span> <span class="o">*</span><span class="n">shape</span><span class="p">,</span> <span class="nb">int</span> <span class="n">nrows</span><span class="p">,</span> <span class="n">nobs</span> <span class="o">=</span> <span class="bp">None</span><span class="p">):</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">==</span> <span class="n">nrows</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Invalid shape for </span><span class="si">%s</span><span class="s"> vector: requires </span><span class="si">%d</span><span class="s"> rows,&#39;</span>
<span class="s">&#39; got </span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">name</span><span class="p">,</span> <span class="n">nrows</span><span class="p">,</span> <span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]))</span>
<span class="k">if</span> <span class="n">nobs</span> <span class="ow">is</span> <span class="ow">not</span> <span class="bp">None</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="ow">in</span> <span class="p">[</span><span class="mf">1</span><span class="p">,</span> <span class="n">nobs</span><span class="p">]:</span>
<span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Invalid time-varying dimension for </span><span class="si">%s</span><span class="s"> vector:&#39;</span>
<span class="s">&#39; requires 1 or </span><span class="si">%d</span><span class="s"> got </span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">name</span><span class="p">,</span> <span class="n">nobs</span><span class="p">,</span> <span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]))</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-14'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-14'>#</a>
</div>
<h1>State Space Representation</h1>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="k">class</span> <span class="nf">dStatespace</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-15'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-15'>#</a>
</div>
<p>dStatespace(obs, design, obs_intercept, obs_cov, transition, state_intercept, selection, state_cov)</p>
<p><em>See Durbin and Koopman (2012), Chapter 4 for all notation</em></p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-16'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-16'>#</a>
</div>
<h3>State space representation</h3>
<p>$$
\begin{align}
y_t &amp; = Z_t \alpha_t + d_t + \varepsilon_t \hspace{3em} &amp; \varepsilon_t &amp; \sim N(0, H_t) \\
\alpha_{t+1} &amp; = T_t \alpha_t + c_t + R_t \eta_t &amp; \eta_t &amp; \sim N(0, Q_t) \\
&amp; &amp; \alpha_1 &amp; \sim N(a_1, P_1)
\end{align}
$$</p>
<p>$y_t$ is $p \times 1$<br />
$\varepsilon_t$ is $p \times 1$<br />
$\alpha_t$ is $m \times 1$<br />
$\eta_t$ is $r \times 1$<br />
$t = 1, \dots, T$</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-17'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-17'>#</a>
</div>
<p><code>nobs</code> $\equiv T$ is the length of the time-series<br />
<code>nendog</code> $\equiv p$ is dimension of observation space<br />
<code>nstates</code> $\equiv m$ is the dimension of the state space<br />
<code>nposdef</code> $\equiv r$ is the dimension of the state shocks<br />
<em>Old notation: T, n, k, g</em></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">nobs</span><span class="p">,</span> <span class="nf">nendog</span><span class="p">,</span> <span class="nf">nstates</span><span class="p">,</span> <span class="nf">nposdef</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-18'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-18'>#</a>
</div>
<p><code>obs</code> $\equiv y_t$ is the <strong>observation vector</strong> $(p \times T)$<br />
<code>design</code> $\equiv Z_t$ is the <strong>design vector</strong> $(p \times m \times T)$<br />
<code>obs_intercept</code> $\equiv d_t$ is the <strong>observation intercept</strong> $(p \times T)$<br />
<code>obs_cov</code> $\equiv H_t$ is the <strong>observation covariance matrix</strong> $(p \times p \times T)$<br />
<code>transition</code> $\equiv T_t$ is the <strong>transition matrix</strong> $(m \times m \times T)$<br />
<code>state_intercept</code> $\equiv c_t$ is the <strong>state intercept</strong> $(m \times T)$<br />
<code>selection</code> $\equiv R_t$ is the <strong>selection matrix</strong> $(m \times r \times T)$<br />
<code>state_cov</code> $\equiv Q_t$ is the <strong>state covariance matrix</strong> $(r \times r \times T)$<br />
<code>selected_state_cov</code> $\equiv R Q_t R'$ is the <strong>selected state covariance matrix</strong> $(m \times m \times T)$<br />
<code>initial_state</code> $\equiv a_1$ is the <strong>initial state mean</strong> $(m \times 1)$<br />
<code>initial_state_cov</code> $\equiv P_1$ is the <strong>initial state covariance matrix</strong> $(m \times m)$</p>
<p>With the exception of <code>obs</code>, these are <em>optionally</em> time-varying. If they are instead time-invariant,
then the dimension of length $T$ is instead of length $1$.</p>
<p><em>Note</em>: the initial vectors' notation 1-indexed as in Durbin and Koopman,
but in the recursions below it will be 0-indexed in the Python arrays.</p>
<p><em>Old notation: y, -, mu, beta_tt_init, P_tt_init</em></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">obs</span><span class="p">,</span> <span class="n">obs_intercept</span><span class="p">,</span> <span class="n">state_intercept</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">initial_state</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">initial_state_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-19'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-19'>#</a>
</div>
<p><em>Old notation: H, R, F, G, Q</em>, G Q<em> G'</em></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">design</span><span class="p">,</span> <span class="n">obs_cov</span><span class="p">,</span> <span class="n">transition</span><span class="p">,</span> <span class="n">selection</span><span class="p">,</span> <span class="n">state_cov</span><span class="p">,</span> <span class="n">selected_state_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-20'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-20'>#</a>
</div>
<p><code>missing</code> is a $(p \times T)$ boolean matrix where a row is a $(p \times 1)$ vector
in which the $i$th position is $1$ if $y_{i,t}$ is to be considered a missing value.<br />
<em>Note:</em> This is created as the output of np.isnan(obs).</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">missing</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-21'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-21'>#</a>
</div>
<p><code>nmissing</code> is an <code>T \times 0</code> integer vector holding the number of <em>missing</em> observations
$p - p_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> [<span class="p">:]</span> <span class="n">nmissing</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-22'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-22'>#</a>
</div>
<p>Flag for a time-invariant model, which requires that <em>all</em> of the
possibly time-varying arrays are time-invariant.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">time_invariant</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-23'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-23'>#</a>
</div>
<p>Flag for initialization.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">initialized</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-24'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-24'>#</a>
</div>
<h3>Kalman filter properties</h3>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-25'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-25'>#</a>
</div>
<p><code>loglikelihood</code> $\equiv \log p(y_t | Y_{t-1})$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">loglikelihood</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-26'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-26'>#</a>
</div>
<p><code>filtered_state</code> $\equiv a_{t|t} = E(\alpha_t | Y_t)$ is the <strong>filtered estimator</strong> of the state $(m \times T)$<br />
<code>predicted_state</code> $\equiv a_{t+1} = E(\alpha_{t+1} | Y_t)$ is the <strong>one-step ahead predictor</strong> of the state $(m \times T-1)$<br />
<code>forecast</code> $\equiv E(y_t|Y_{t-1})$ is the <strong>forecast</strong> of the next observation $(p \times T)$ <br />
<code>forecast_error</code> $\equiv v_t = y_t - E(y_t|Y_{t-1})$ is the <strong>one-step ahead forecast error</strong> of the next observation $(p \times T)$ </p>
<p><em>Note</em>: Actual values in <code>filtered_state</code> will be from 1 to <code>nobs</code>+1. Actual
values in <code>predicted_state</code> will be from 0 to <code>nobs</code>+1 because the initialization
is copied over to the zeroth entry, and similar for the covariances, below.</p>
<p><em>Old notation: beta_tt, beta_tt1, y_tt1, eta_tt1</em></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">filtered_state</span><span class="p">,</span> <span class="n">predicted_state</span><span class="p">,</span> <span class="n">forecast</span><span class="p">,</span> <span class="n">forecast_error</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-27'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-27'>#</a>
</div>
<p><code>filtered_state_cov</code> $\equiv P_{t|t} = Var(\alpha_t | Y_t)$ is the <strong>filtered state covariance matrix</strong> $(m \times m \times T)$<br />
<code>predicted_state_cov</code> $\equiv P_{t+1} = Var(\alpha_{t+1} | Y_t)$ is the <strong>predicted state covariance matrix</strong> $(m \times m \times T)$<br />
<code>forecast_error_cov</code> $\equiv F_t = Var(v_t | Y_{t-1})$ is the <strong>forecast error covariance matrix</strong> $(p \times p \times T)$ </p>
<p><em>Old notation: P_tt, P_tt1, f_tt1</em></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">filtered_state_cov</span><span class="p">,</span> <span class="n">predicted_state_cov</span><span class="p">,</span> <span class="n">forecast_error_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-28'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-28'>#</a>
</div>
<h3>Initialize state space model</h3>
<p><em>Note</em>: The initial state and state covariance matrix must be provided.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">obs</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">design</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">obs_intercept</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">obs_cov</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">transition</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">state_intercept</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">selection</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:,:]</span> <span class="n">state_cov</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-29'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-29'>#</a>
</div>
<p>Local variables</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim1</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">3</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-30'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-30'>#</a>
</div>
<h4>State space representation variables</h4>
<p><strong>Note</strong>: these arrays share data with the versions defined in
Python and passed to this constructor, so if they are updated in
Python they will also be updated here.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">obs</span> <span class="o">=</span> <span class="n">obs</span>
<span class="bp">self</span><span class="o">.</span><span class="n">design</span> <span class="o">=</span> <span class="n">design</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_intercept</span> <span class="o">=</span> <span class="n">obs_intercept</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_cov</span> <span class="o">=</span> <span class="n">obs_cov</span>
<span class="bp">self</span><span class="o">.</span><span class="n">transition</span> <span class="o">=</span> <span class="n">transition</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_intercept</span> <span class="o">=</span> <span class="n">state_intercept</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selection</span> <span class="o">=</span> <span class="n">selection</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span> <span class="o">=</span> <span class="n">state_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-31'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-31'>#</a>
</div>
<p>Dimensions</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">=</span> <span class="n">obs</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nstates</span> <span class="o">=</span> <span class="n">selection</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span> <span class="o">=</span> <span class="n">selection</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nobs</span> <span class="o">=</span> <span class="n">obs</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-32'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-32'>#</a>
</div>
<h4>Validate matrix dimensions</h4>
<p>Make sure that the given state-space matrices have consistent sizes</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">validate_matrix_shape</span><span class="p">(</span><span class="s">&#39;design&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">design</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span>
<span class="n">validate_vector_shape</span><span class="p">(</span><span class="s">&#39;observation intercept&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">obs_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span>
<span class="n">validate_matrix_shape</span><span class="p">(</span><span class="s">&#39;observation covariance matrix&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">obs_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span>
<span class="n">validate_matrix_shape</span><span class="p">(</span><span class="s">&#39;transition&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">transition</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span>
<span class="n">validate_vector_shape</span><span class="p">(</span><span class="s">&#39;state intercept&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">state_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span>
<span class="n">validate_matrix_shape</span><span class="p">(</span><span class="s">&#39;state covariance matrix&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-33'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-33'>#</a>
</div>
<p>Check for a time-invariant model</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">time_invariant</span> <span class="o">=</span> <span class="p">(</span>
<span class="bp">self</span><span class="o">.</span><span class="n">design</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">transition</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selection</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span> <span class="ow">and</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">==</span> <span class="mf">1</span>
<span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-34'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-34'>#</a>
</div>
<p>Set the flag for initialization to be false</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">initialized</span> <span class="o">=</span> <span class="bp">False</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-35'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-35'>#</a>
</div>
<h4>Allocate arrays for calculations</h4>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-36'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-36'>#</a>
</div>
<p>Selected state covariance matrix</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim3</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1</span><span class="p">;</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-37'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-37'>#</a>
</div>
<p>(we only allocate memory for time-varying array if necessary)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span> <span class="ow">or</span> <span class="bp">self</span><span class="o">.</span><span class="n">selection</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span>
<span class="n">dim3</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_state_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">3</span><span class="p">,</span> <span class="n">dim3</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-38'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-38'>#</a>
</div>
<p>Arrays for Kalman filter output</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">filtered_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="o">+</span><span class="mf">1</span>
<span class="bp">self</span><span class="o">.</span><span class="n">predicted_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim3</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">filtered_state_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">3</span><span class="p">,</span> <span class="n">dim3</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim3</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="o">+</span><span class="mf">1</span>
<span class="bp">self</span><span class="o">.</span><span class="n">predicted_state_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">3</span><span class="p">,</span> <span class="n">dim3</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim3</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">3</span><span class="p">,</span> <span class="n">dim3</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nobs</span>
<span class="bp">self</span><span class="o">.</span><span class="n">loglikelihood</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-39'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-39'>#</a>
</div>
<p>Handle missing data</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">missing</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">obs</span><span class="p">),</span> <span class="n">dtype</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">int32</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="s">&quot;F&quot;</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nmissing</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">missing</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mf">0</span><span class="p">),</span> <span class="n">dtype</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">int32</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-40'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-40'>#</a>
</div>
<h2>Initialize: known values</h2>
<p>Initialize the filter with specific values, assumed to be known with
certainty or else as filled with parameters from a maximum likelihood
estimation run.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">initialize_known</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[:]</span> <span class="n">initial_state</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="p">[::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">initial_state_cov</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-41'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-41'>#</a>
</div>
<p>initialize_known(initial_state, initial_state_cov)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">validate_vector_shape</span><span class="p">(</span><span class="s">&#39;inital state&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">initial_state</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">None</span><span class="p">)</span>
<span class="n">validate_matrix_shape</span><span class="p">(</span><span class="s">&#39;initial state covariance&#39;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">initial_state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">0</span><span class="p">],</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">None</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state</span> <span class="o">=</span> <span class="n">initial_state</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state_cov</span> <span class="o">=</span> <span class="n">initial_state_cov</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initialized</span> <span class="o">=</span> <span class="bp">True</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-42'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-42'>#</a>
</div>
<h2>Initialize: approximate diffuse priors</h2>
<p>Durbin and Koopman note that this initialization should only be coupled
with the standard Kalman filter for "approximate exploratory work" and
can lead to "large rounding errors" (p. 125).</p>
<p><em>Note:</em> see Durbin and Koopman section 5.6.1</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">initialize_approximate_diffuse</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">variance</span><span class="o">=</span><span class="mf">1e2</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-43'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-43'>#</a>
</div>
<p>initialize_approximate_diffuse(variance=1e2)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">npy_intp</span> <span class="kt">dim</span>[1]
<span class="n">dim</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span><span class="o">.</span><span class="n">T</span> <span class="o">*</span> <span class="n">variance</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initialized</span> <span class="o">=</span> <span class="bp">True</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-44'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-44'>#</a>
</div>
<h2>Initialize: stationary process</h2>
<p><em>Note:</em> see Durbin and Koopman section 5.6.2</p>
<p>TODO improve efficiency with direct BLAS / LAPACK calls</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">initialize_stationary</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-45'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-45'>#</a>
</div>
<p>initialize_stationary()</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">npy_intp</span> <span class="kt">dim</span>[1]
<span class="k">from</span> <span class="nn">scipy.linalg</span> <span class="k">import</span> <span class="n">solve_discrete_lyapunov</span>
<span class="n">dim</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state_cov</span> <span class="o">=</span> <span class="n">solve_discrete_lyapunov</span><span class="p">(</span>
<span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">transition</span><span class="p">[:,:,</span><span class="mf">0</span><span class="p">],</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">),</span>
<span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span><span class="p">[:,:,</span><span class="mf">0</span><span class="p">],</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">float</span><span class="p">)</span>
<span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initialized</span> <span class="o">=</span> <span class="bp">True</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-46'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-46'>#</a>
</div>
<h2>Kalman filter Routines</h2>
<p>The following functions are the workhorse functions for the Kalman filter.
They represent four distinct but very general phases of the Kalman filtering
operations.</p>
<p>Their argument is an object of class ?KalmanFilter, which is a stateful
representation of the recursive filter. For this reason, the below functions
work almost exclusively through <em>side-effects</em> and most return void.
See the Kalman filter class documentation for further discussion.</p>
<p>They are defined this way so that the actual filtering process can select
whichever filter type is appropriate for the given time period. For example,
in the case of state space models with non-stationary components, the filter
should begin with the exact initial Kalman filter routines but after some
number of time periods will transition to the conventional Kalman filter
routines.</p>
<p>Below, <code>&lt;filter type&gt;</code> will refer to one of the following:</p>
<ul>
<li><code>conventional</code> - the conventional Kalman filter</li>
</ul>
<p>Other filter types (e.g. <code>exact_initial</code>, <code>augmented</code>, etc.) may be added in
the future.</p>
<p><code>forecast_&lt;filter type&gt;</code> generates the forecast, forecast error $v_t$ and
forecast error covariance matrix $F_t$<br />
<code>updating_&lt;filter type&gt;</code> is the updating step of the Kalman filter, and
generates the filtered state $a_{t|t}$ and covariance matrix $P_{t|t}$<br />
<code>prediction_&lt;filter type&gt;</code> is the prediction step of the Kalman filter, and
generates the predicted state $a_{t+1}$ and covariance matrix $P_{t+1}$.
<code>loglikelihood_&lt;filter type&gt;</code> calculates the loglikelihood for $y_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-47'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-47'>#</a>
</div>
<h3>Conventional Kalman filter</h3>
<p>The following are the above routines as defined in the conventional Kalman
filter.</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="kt">void</span> <span class="nf">dforecast_conventional</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-48'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-48'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">gamma</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.0</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-49'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-49'>#</a>
</div>
<h4>Forecast for time t</h4>
<p><code>forecast</code> $= Z_t a_t + d_t$</p>
<p><em>Note</em>: $a_t$ is given from the initialization (for $t = 0$) or
from the previous iteration of the filter (for $t &gt; 0$).</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-50'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-50'>#</a>
</div>
<p>$\# = d_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">obs_intercept</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-51'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-51'>#</a>
</div>
<p><code>forecast</code> $= 1.0 * Z_t a_t + 1.0 * \#$<br />
$(p \times 1) = (p \times m) (m \times 1) + (p \times 1)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemv</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">input_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-52'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-52'>#</a>
</div>
<h4>Forecast error for time t</h4>
<p><code>forecast_error</code> $\equiv v_t = y_t -$ <code>forecast</code></p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-53'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-53'>#</a>
</div>
<p>$\# = y_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">obs</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-54'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-54'>#</a>
</div>
<p>$v_t = -1.0 * $ <code>forecast</code> $ + \#$
$(p \times 1) = (p \times 1) + (p \times 1)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">daxpy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">gamma</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-55'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-55'>#</a>
</div>
<p><em>Intermediate calculation</em> (used just below and then once more)<br />
<code>tmp1</code> array used here, dimension $(m \times p)$<br />
$\#_1 = P_t Z_t'$<br />
$(m \times p) = (m \times m) (p \times m)'$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;T&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">input_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp1</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-56'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-56'>#</a>
</div>
<h4>Forecast error covariance matrix for time t</h4>
<p>$F_t \equiv Z_t P_t Z_t' + H_t$</p>
<p><em>Note</em>: this and does nothing at all to <code>forecast_error_cov</code> if
converged == True</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-57'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-57'>#</a>
</div>
<p>$\# = H_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog2</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">obs_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-58'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-58'>#</a>
</div>
<p>$F_t = 1.0 * Z_t \#_1 + 1.0 * \#$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">tmp1</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">dupdating_conventional</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-59'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-59'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">gamma</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.0</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-60'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-60'>#</a>
</div>
<h4>Filtered state for time t</h4>
<p>$a_{t|t} = a_t + P_t Z_t' F_t^{-1} v_t$<br />
$a_{t|t} = 1.0 * \#_1 \#_2 + 1.0 a_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">input_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dgemv</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp1</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-61'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-61'>#</a>
</div>
<h4>Filtered state covariance for time t</h4>
<p>$P_{t|t} = P_t - P_t Z_t' F_t^{-1} Z_t P_t$<br />
$P_{t|t} = P_t - \#_1 \#_3 P_t$ </p>
<p><em>Note</em>: this and does nothing at all to <code>filtered_state_cov</code> if
converged == True</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span>
<span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">input_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-62'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-62'>#</a>
</div>
<p><code>tmp0</code> array used here, dimension $(m \times m)$<br />
$\#_0 = 1.0 * \#_1 \#_3$<br />
$(m \times m) = (m \times p) (p \times m)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp1</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-63'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-63'>#</a>
</div>
<p>$P_{t|t} = - 1.0 * \# P_t + 1.0 * P_t$<br />
$(m \times m) = (m \times m) (m \times m) + (m \times m)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">gamma</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">input_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">dprediction_conventional</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-64'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-64'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">gamma</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.0</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-65'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-65'>#</a>
</div>
<h4>Predicted state for time t+1</h4>
<p>$a_{t+1} = T_t a_{t|t} + c_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">state_intercept</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">predicted_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dgemv</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">transition</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">predicted_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-66'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-66'>#</a>
</div>
<h4>Predicted state covariance matrix for time t+1</h4>
<p>$P_{t+1} = T_t P_{t|t} T_t' + Q_t^*$</p>
<p><em>Note</em>: this and does nothing at all to <code>predicted_state_cov</code> if
converged == True</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span>
<span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">selected_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-67'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-67'>#</a>
</div>
<p><code>tmp0</code> array used here, dimension $(m \times m)$ </p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-68'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-68'>#</a>
</div>
<p>$\#<em>0 = T_t P</em>{t|t} $</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-69'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-69'>#</a>
</div>
<p>$(m \times m) = (m \times m) (m \times m)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">transition</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-70'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-70'>#</a>
</div>
<p>$P_{t+1} = 1.0 \#_0 T_t' + 1.0 \#$<br />
$(m \times m) = (m \times m) (m \times m) + (m \times m)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;T&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">transition</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dloglikelihood_conventional</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-71'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-71'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">loglikelihood</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">loglikelihood</span> <span class="o">=</span> <span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="p">(</span><span class="n">dlog</span><span class="p">(</span><span class="mf">2</span><span class="o">*</span><span class="n">NPY_PI</span><span class="o">*</span><span class="n">determinant</span><span class="p">))</span>
<span class="n">loglikelihood</span> <span class="o">=</span> <span class="n">loglikelihood</span> <span class="o">-</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">ddot</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="k">return</span> <span class="n">loglikelihood</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-72'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-72'>#</a>
</div>
<h2>Forecast error covariance inversion</h2>
<p>The following are routines that can calculate the inverse of the forecast
error covariance matrix (defined in <code>forecast_&lt;filter type&gt;</code>).</p>
<p>These routines are aware of the possibility that the Kalman filter may have
converged to a steady state, in which case they do not need to perform the
inversion or calculate the determinant.</p>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dinverse_univariate</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-73'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-73'>#</a>
</div>
<p>Factorize the forecast error covariance matrix using simple division
in the case that the observations are univariate.</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
<h4>Intermediate values</h4>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">scalar</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-74'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-74'>#</a>
</div>
<p>Take the inverse of the forecast error covariance matrix</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span>
<span class="n">determinant</span> <span class="o">=</span> <span class="n">dabs</span><span class="p">(</span><span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">])</span>
<span class="n">scalar</span> <span class="o">=</span> <span class="mf">1.0</span> <span class="o">/</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">scalar</span> <span class="o">*</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendogstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dscal</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendogstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">scalar</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dfactorize_cholesky</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-75'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-75'>#</a>
</div>
<p>Factorize the forecast error covariance matrix using a Cholesky
decomposition. Called by either of the <code>solve_cholesky</code> or
<code>invert_cholesky</code> routines.</p>
<p>Requires a positive definite matrix, but is faster than an LU
decomposition.</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">i</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span>
<span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog2</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dpotrf</span><span class="p">(</span><span class="s">&quot;U&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span>
<span class="k">if</span> <span class="n">info</span> <span class="o">&gt;</span> <span class="mf">0</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-76'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-76'>#</a>
</div>
<p>TODO catch and re-raise this exception in main loop so that e.g.
the time subscript can be provided as well.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">warnings</span><span class="o">.</span><span class="n">warn</span><span class="p">(</span><span class="s">&#39;Non-positive-definite forecast error&#39;</span>
<span class="s">&#39; covariance matrix encountered.&#39;</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-77'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-77'>#</a>
</div>
<p>Calculate the determinant (just the squared product of the
diagonals, in the Cholesky decomposition case)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">):</span>
<span class="n">determinant</span> <span class="o">=</span> <span class="n">determinant</span> <span class="o">*</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">i</span><span class="p">]</span>
<span class="n">determinant</span> <span class="o">=</span> <span class="n">determinant</span><span class="o">**</span><span class="mf">2</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dfactorize_lu</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-78'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-78'>#</a>
</div>
<p>Factorize the forecast error covariance matrix using an LU
decomposition. Called by either of the <code>solve_lu</code> or <code>invert_lu</code>
routines.</p>
<p>Is slower than a Cholesky decomposition, but does not require a
positive definite matrix.</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">i</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-79'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-79'>#</a>
</div>
<p>Perform LU decomposition into <code>forecast_error_fac</code></p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog2</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dgetrf</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_ipiv</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-80'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-80'>#</a>
</div>
<p>Calculate the determinant (product of the diagonals, but with
sign modifications according to the permutation matrix) </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="mf">1</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">):</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_ipiv</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="n">i</span><span class="o">+</span><span class="mf">1</span><span class="p">:</span>
<span class="n">determinant</span> <span class="o">*=</span> <span class="o">-</span><span class="mf">1</span><span class="o">*</span><span class="n">kfilter</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">i</span><span class="p">]</span>
<span class="k">else</span><span class="p">:</span>
<span class="n">determinant</span> <span class="o">*=</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">i</span><span class="p">]</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dinverse_cholesky</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-81'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-81'>#</a>
</div>
<p>inverse_cholesky(self, determinant)</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="nb">int</span> <span class="n">i</span><span class="p">,</span> <span class="n">j</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-82'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-82'>#</a>
</div>
<p>Perform the Cholesky decomposition and get the determinant</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="n">dfactorize_cholesky</span><span class="p">(</span><span class="n">kfilter</span><span class="p">,</span> <span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-83'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-83'>#</a>
</div>
<p>Continue taking the inverse</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dpotri</span><span class="p">(</span><span class="s">&quot;U&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-84'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-84'>#</a>
</div>
<p>?potri only fills in the upper triangle of the symmetric array, and
since the ?symm and ?symv routines are not available as of scipy
0.11.0, we can't use them, so we must fill in the lower triangle
by hand</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">):</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">i</span><span class="p">):</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-85'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-85'>#</a>
</div>
<p>Get <code>tmp2</code> and <code>tmp3</code> via matrix multiplications</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-86'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-86'>#</a>
</div>
<p><code>tmp2</code> array used here, dimension $(p \times 1)$<br />
$\#_2 = F_t^{-1} v_t$<br />
dsymv("U", &amp;kfilter.nendog, &amp;alpha, kfilter.forecast_error_fac, &amp;kfilter.nendog,
kfilter.forecast_error, &amp;inc, &amp;beta, kfilter.tmp2, &amp;inc)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemv</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-87'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-87'>#</a>
</div>
<p><code>tmp3</code> array used here, dimension $(p \times m)$<br />
$\#_3 = F_t^{-1} Z_t$
dsymm("L", "U", &amp;kfilter.nendog, &amp;kfilter.nstates,
&amp;alpha, kfilter.forecast_error_fac, &amp;kfilter.nendog,
kfilter.design, &amp;kfilter.nendog,
&amp;beta, kfilter.tmp3, &amp;kfilter.nendog)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">)</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dinverse_lu</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-88'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-88'>#</a>
</div>
<p>inverse_cholesky(self, determinant)</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-89'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-89'>#</a>
</div>
<p>Perform the Cholesky decomposition and get the determinant</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="n">dfactorize_lu</span><span class="p">(</span><span class="n">kfilter</span><span class="p">,</span> <span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-90'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-90'>#</a>
</div>
<p>Continue taking the inverse</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgetri</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_ipiv</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_work</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">ldwork</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-91'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-91'>#</a>
</div>
<p>Get <code>tmp2</code> and <code>tmp3</code> via matrix multiplications</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-92'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-92'>#</a>
</div>
<p><code>tmp2</code> array used here, dimension $(p \times 1)$<br />
$\#_2 = F_t^{-1} v_t$ </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemv</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-93'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-93'>#</a>
</div>
<p><code>tmp3</code> array used here, dimension $(p \times m)$<br />
$\#_3 = F_t^{-1} Z_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">)</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dsolve_cholesky</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-94'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-94'>#</a>
</div>
<p>solve_cholesky(self, determinant)</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-95'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-95'>#</a>
</div>
<p>Perform the Cholesky decomposition and get the determinant</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="n">dfactorize_cholesky</span><span class="p">(</span><span class="n">kfilter</span><span class="p">,</span> <span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-96'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-96'>#</a>
</div>
<p>Solve the linear systems<br />
<code>tmp2</code> array used here, dimension $(p \times 1)$<br />
$F_t \#_2 = v_t$ </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dpotrs</span><span class="p">(</span><span class="s">&quot;U&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-97'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-97'>#</a>
</div>
<p><code>tmp3</code> array used here, dimension $(p \times m)$<br />
$F_t \#_3 = Z_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendogstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dpotrs</span><span class="p">(</span><span class="s">&quot;U&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span>
<span class="k">return</span> <span class="n">determinant</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">dsolve_lu</span><span class="p">(</span><span class="n">dKalmanFilter</span> <span class="n">kfilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">determinant</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-98'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-98'>#</a>
</div>
<p>inverse_cholesky(self, determinant)</p>
<p>If the model has converged to a steady-state, this is a NOOP and simply
returns the determinant that was passed in.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">info</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-99'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-99'>#</a>
</div>
<p>Perform the Cholesky decomposition and get the determinant</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">determinant</span> <span class="o">=</span> <span class="n">dfactorize_lu</span><span class="p">(</span><span class="n">kfilter</span><span class="p">,</span> <span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-100'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-100'>#</a>
</div>
<p>Solve the linear systems<br />
<code>tmp2</code> array used here, dimension $(p \times 1)$<br />
$F_t \#_2 = v_t$ </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dgetrs</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_ipiv</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-101'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-101'>#</a>
</div>
<p><code>tmp3</code> array used here, dimension $(p \times m)$<br />
$F_t \#_3 = Z_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendogstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">design</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dgetrs</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_fac</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="n">kfilter</span><span class="o">.</span><span class="n">forecast_error_ipiv</span><span class="p">,</span> <span class="n">kfilter</span><span class="o">.</span><span class="n">tmp3</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">kfilter</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">info</span><span class="p">)</span>
<span class="k">return</span> <span class="n">determinant</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-102'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-102'>#</a>
</div>
<h2>Kalman filter</h2>
</div>
<div class='code'>
<div class="highlight"><pre><span class="k">cdef</span> <span class="k">class</span> <span class="nf">dKalmanFilter</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-103'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-103'>#</a>
</div>
<p>dKalmanFilter(model, filter=FILTER_CONVENTIONAL, inversion_method=INVERT_UNIVARIATE | SOLVE_CHOLESKY, stability_method=STABILITY_FORCE_SYMMETRY, tolerance=1e-19)</p>
<p>A representation of the Kalman filter recursions.</p>
<p>While the filter is mathematically represented as a recursion, it is here
translated into Python as a stateful iterator.</p>
<p>Because there are actually several types of Kalman filter depending on the
state space model of interest, this class only handles the <em>iteration</em>
aspect of filtering, and delegates the actual operations to four general
workhorse routines, which can be implemented separately for each type of
Kalman filter.</p>
<p>In order to maintain a consistent interface, and because these four general
routines may be quite different across filter types, their argument is only
the stateful ?KalmanFilter object. Furthermore, in order to allow the
different types of filter to substitute alternate matrices, this class
defines a set of pointers to the various state space arrays and the
filtering output arrays.</p>
<p>For example, handling missing observations requires not only substituting
<code>obs</code>, <code>design</code>, and <code>obs_cov</code> matrices, but the new matrices actually have
different dimensions than the originals. This can be flexibly accomodated
simply by replacing e.g. the <code>obs</code> pointer to the substituted <code>obs</code> array
and replacing <code>nendog</code> for that iteration. Then in the next iteration, when
the <code>obs</code> vector may be missing different elements (or none at all), it can
again be redefined.</p>
<p>Each iteration of the filter (see <code>__next__</code>) proceeds in a number of
steps.</p>
<p><code>initialize_object_pointers</code> initializes pointers to current-iteration
objects (i.e. the state space arrays and filter output arrays). </p>
<p><code>initialize_function_pointers</code> initializes pointers to the appropriate
Kalman filtering routines (i.e. <code>forecast_conventional</code> or
<code>forecast_exact_initial</code>, etc.). </p>
<p><code>select_arrays</code> converts the base arrays into "selected" arrays using
selection matrices. In particular, it handles the state covariance matrix
and redefined matrices based on missing values. </p>
<p><code>post_convergence</code> handles copying arrays from time $t-1$ to time $t$ when
the Kalman filter has converged and they don't need to be re-calculated. </p>
<p><code>forecasting</code> calls the Kalman filter <code>forcasting_&lt;filter type&gt;</code> routine</p>
<p><code>inversion</code> calls the appropriate function to invert the forecast error
covariance matrix. </p>
<p><code>updating</code> calls the Kalman filter <code>updating_&lt;filter type&gt;</code> routine</p>
<p><code>loglikelihood</code> calls the Kalman filter <code>loglikelihood_&lt;filter type&gt;</code> routine</p>
<p><code>prediction</code> calls the Kalman filter <code>predcition_&lt;filter type&gt;</code> routine</p>
<p><code>numerical_stability</code> performs end-of-iteration tasks to improve the numerical
stability of the filter </p>
<p><code>check_convergence</code> checks for convergence of the filter to steady-state.</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-104'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-104'>#</a>
</div>
<h3>Statespace model</h3>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">dStatespace</span> <span class="nf">model</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-105'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-105'>#</a>
</div>
<h3>Filter parameters</h3>
<p>Holds the time-iteration state of the filter<br />
<em>Note</em>: must be changed using the <code>seek</code> method</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">t</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-106'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-106'>#</a>
</div>
<p>Holds the tolerance parameter for convergence</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">public</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">tolerance</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-107'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-107'>#</a>
</div>
<p>Holds the convergence to steady-state status of the filter
<em>Note</em>: is by default reset each time <code>seek</code> is called</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">converged</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-108'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-108'>#</a>
</div>
<p>Holds whether or not the model is time-invariant
<em>Note</em>: is by default reset each time <code>seek</code> is called</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> <span class="nf">time_invariant</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-109'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-109'>#</a>
</div>
<p>The Kalman filter procedure to use </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">public</span> <span class="kt">int</span> <span class="nf">filter_method</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-110'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-110'>#</a>
</div>
<p>The method by which the terms using the inverse of the forecast
error covariance matrix are solved.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">public</span> <span class="kt">int</span> <span class="nf">inversion_method</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-111'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-111'>#</a>
</div>
<p>Methods to improve numerical stability</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">public</span> <span class="kt">int</span> <span class="nf">stability_method</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-112'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-112'>#</a>
</div>
<h3>Temporary arrays</h3>
<p>These matrices are used to temporarily hold selected observation vectors,
design matrices, and observation covariance matrices in the case of
missing data. </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">selected_obs</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">_selected_design</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">_selected_obs_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-113'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-113'>#</a>
</div>
<p>The following are contiguous memory segments which are then used to
store the data in the above matrices.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">selected_design</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">selected_obs_cov</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-114'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-114'>#</a>
</div>
<p><code>forecast_error_fac</code> is a forecast error covariance matrix <strong>factorization</strong> $(p \times p)$.
Depending on the method for handling the inverse of the forecast error covariance matrix, it may be:
- a Cholesky factorization if <code>cholesky_solve</code> is used
- an inverse calculated via Cholesky factorization if <code>cholesky_inverse</code> is used
- an LU factorization if <code>lu_solve</code> is used
- an inverse calculated via LU factorization if <code>lu_inverse</code> is used</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">_forecast_error_fac</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-115'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-115'>#</a>
</div>
<p><code>forecast_error_ipiv</code> holds pivot indices if an LU decomposition is used</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">int</span> [<span class="p">:]</span> <span class="n">_forecast_error_ipiv</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-116'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-116'>#</a>
</div>
<p><code>forecast_error_work</code> is a work array for matrix inversion if an LU
decomposition is used</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">_forecast_error_work</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-117'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-117'>#</a>
</div>
<p>These hold the memory allocations of the unnamed temporary arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">::</span><span class="mf">1</span><span class="p">,:]</span> <span class="n">_tmp0</span><span class="p">,</span> <span class="n">_tmp1</span><span class="p">,</span> <span class="n">_tmp3</span>
<span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> [<span class="p">:]</span> <span class="n">_tmp2</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-118'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-118'>#</a>
</div>
<p>Holds the determinant across calculations (this is done because after
convergence, it doesn't need to be re-calculated anymore)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kr">readonly</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">determinant</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-119'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-119'>#</a>
</div>
<h3>Pointers to current-iteration arrays</h3>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">obs</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">design</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">obs_intercept</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">obs_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">transition</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">state_intercept</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">selection</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">selected_state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">initial_state</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">initial_state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">input_state</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">input_state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">forecast</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">forecast_error</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">forecast_error_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">filtered_state</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">filtered_state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">predicted_state</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">predicted_state_cov</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">forecast_error_fac</span>
<span class="k">cdef</span> <span class="kt">int</span> * <span class="nf">forecast_error_ipiv</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">forecast_error_work</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">tmp0</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">tmp1</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">tmp2</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> * <span class="nf">tmp3</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-120'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-120'>#</a>
</div>
<h3>Pointers to current-iteration Kalman filtering functions</h3>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="nf">void</span> <span class="p">(</span><span class="o">*</span><span class="n">forecasting</span><span class="p">)(</span>
<span class="n">dKalmanFilter</span>
<span class="p">)</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="nf">float64_t</span> <span class="p">(</span><span class="o">*</span><span class="n">inversion</span><span class="p">)(</span>
<span class="n">dKalmanFilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span>
<span class="p">)</span>
<span class="k">cdef</span> <span class="nf">void</span> <span class="p">(</span><span class="o">*</span><span class="n">updating</span><span class="p">)(</span>
<span class="n">dKalmanFilter</span>
<span class="p">)</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="nf">float64_t</span> <span class="p">(</span><span class="o">*</span><span class="n">loglikelihood</span><span class="p">)(</span>
<span class="n">dKalmanFilter</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">float64_t</span>
<span class="p">)</span>
<span class="k">cdef</span> <span class="nf">void</span> <span class="p">(</span><span class="o">*</span><span class="n">prediction</span><span class="p">)(</span>
<span class="n">dKalmanFilter</span>
<span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-121'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-121'>#</a>
</div>
<h3>Define some constants</h3>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kt">int</span> <span class="nf">nendog</span><span class="p">,</span> <span class="nf">nstates</span><span class="p">,</span> <span class="nf">nposdef</span><span class="p">,</span> <span class="nf">nendog2</span><span class="p">,</span> <span class="nf">nstates2</span><span class="p">,</span> <span class="nf">nendogstates</span><span class="p">,</span> <span class="nf">ldwork</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-122'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-122'>#</a>
</div>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span>
<span class="n">dStatespace</span> <span class="n">model</span><span class="p">,</span>
<span class="nb">int</span> <span class="n">filter_method</span><span class="o">=</span><span class="n">FILTER_CONVENTIONAL</span><span class="p">,</span>
<span class="nb">int</span> <span class="n">inversion_method</span><span class="o">=</span><span class="n">INVERT_UNIVARIATE</span> <span class="o">|</span> <span class="n">SOLVE_CHOLESKY</span><span class="p">,</span>
<span class="nb">int</span> <span class="n">stability_method</span><span class="o">=</span><span class="n">STABILITY_FORCE_SYMMETRY</span><span class="p">,</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">tolerance</span><span class="o">=</span><span class="mf">1e-19</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-123'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-123'>#</a>
</div>
<p>Local variables</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim1</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim3</span><span class="p">[</span><span class="mf">3</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-124'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-124'>#</a>
</div>
<p>Save the model</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">model</span> <span class="o">=</span> <span class="n">model</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-125'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-125'>#</a>
</div>
<p>Initialize filter parameters</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">tolerance</span> <span class="o">=</span> <span class="n">tolerance</span>
<span class="k">if</span> <span class="ow">not</span> <span class="n">filter_method</span> <span class="o">==</span> <span class="n">FILTER_CONVENTIONAL</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span><span class="s">&quot;Only the conventional Kalman filter is currently implemented&quot;</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">filter_method</span> <span class="o">=</span> <span class="n">filter_method</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">=</span> <span class="n">inversion_method</span>
<span class="bp">self</span><span class="o">.</span><span class="n">stability_method</span> <span class="o">=</span> <span class="n">stability_method</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-126'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-126'>#</a>
</div>
<p>Initialize the constant values</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">time_invariant</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">time_invariant</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nstates</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nposdef</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog2</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="o">**</span><span class="mf">2</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="o">**</span><span class="mf">2</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendogstates</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-127'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-127'>#</a>
</div>
<p>TODO replace with optimal work array size</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">ldwork</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-128'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-128'>#</a>
</div>
<p>Arrays for temporary calculations
<em>Note</em>: in math notation below, a $\#$ will represent a generic
temporary array, and a $\#_i$ will represent a named temporary array.</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-129'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-129'>#</a>
</div>
<p>Arrays related to matrix factorizations / inverses</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_fac</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_fac</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_fac</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">0</span><span class="p">]</span>
<span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">ldwork</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">ldwork</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_work</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_work</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_work</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">0</span><span class="p">]</span>
<span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_ipiv</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_INT</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_ipiv</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_forecast_error_ipiv</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-130'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-130'>#</a>
</div>
<p>Holds arrays of dimension $(m \times m)$ and $(m \times r)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_tmp0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_tmp0</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-131'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-131'>#</a>
</div>
<p>Holds arrays of dimension $(m \times p)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_tmp1</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">tmp1</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_tmp1</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-132'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-132'>#</a>
</div>
<p>Holds arrays of dimension $(p \times 1)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_tmp2</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">tmp2</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_tmp2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-133'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-133'>#</a>
</div>
<p>Holds arrays of dimension $(p \times m)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim2</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">_tmp3</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">2</span><span class="p">,</span> <span class="n">dim2</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">tmp3</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">_tmp3</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-134'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-134'>#</a>
</div>
<p>Arrays for missing data</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_obs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_design</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span>
<span class="n">dim1</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog2</span><span class="p">;</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_obs_cov</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">PyArray_ZEROS</span><span class="p">(</span><span class="mf">1</span><span class="p">,</span> <span class="n">dim1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">NPY_FLOAT64</span><span class="p">,</span> <span class="n">FORTRAN</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-135'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-135'>#</a>
</div>
<p>Initialize time and convergence status</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">seek</span><span class="p">(</span><span class="mf">0</span><span class="p">)</span>
<span class="k">cpdef</span> <span class="nf">seek</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">unsigned</span> <span class="nb">int</span> <span class="n">t</span><span class="p">,</span> <span class="nb">int</span> <span class="n">reset_convergence</span> <span class="o">=</span> <span class="bp">True</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-136'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-136'>#</a>
</div>
<p>seek(self, t, reset_convergence = True)</p>
<p>Change the time-state of the filter</p>
<p>Is usually called to reset the filter to the beginning.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="n">t</span> <span class="o">&gt;=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nobs</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">IndexError</span><span class="p">(</span><span class="s">&quot;Observation index out of range&quot;</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="n">reset_convergence</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">converged</span> <span class="o">=</span> <span class="mf">0</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-137'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-137'>#</a>
</div>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">__iter__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
<span class="k">return</span> <span class="bp">self</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-138'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-138'>#</a>
</div>
<p>Iterate the filter across the entire set of observations.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">__call__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-139'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-139'>#</a>
</div>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">cdef</span> <span class="kt">int</span> <span class="nf">i</span>
<span class="bp">self</span><span class="o">.</span><span class="n">seek</span><span class="p">(</span><span class="mf">0</span><span class="p">,</span> <span class="bp">True</span><span class="p">)</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nobs</span><span class="p">):</span>
<span class="nb">next</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-140'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-140'>#</a>
</div>
<p>Perform an iteration of the Kalman filter</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">def</span> <span class="nf">__next__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-141'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-141'>#</a>
</div>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-142'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-142'>#</a>
</div>
<p>Get time subscript, and stop the iterator if at the end</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nobs</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">StopIteration</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-143'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-143'>#</a>
</div>
<p>Initialize pointers to current-iteration objects</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">initialize_object_pointers</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-144'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-144'>#</a>
</div>
<p>Initialize pointers to appropriate Kalman filtering functions</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">initialize_function_pointers</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-145'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-145'>#</a>
</div>
<p>Convert base arrays into "selected" arrays<br />
- State covariance matrix? $Q_t \to R_t Q_t R_t`$
- Missing values: $y_t \to W_t y_t$, $Z_t \to W_t Z_t$, $H_t \to W_t H_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">select_arrays</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-146'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-146'>#</a>
</div>
<p>Post-convergence: copy previous iteration arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">post_convergence</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-147'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-147'>#</a>
</div>
<p>Form forecasts</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">forecasting</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-148'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-148'>#</a>
</div>
<p>Perform <code>forecast_error_cov</code> inversion (or decomposition)</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">determinant</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-149'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-149'>#</a>
</div>
<p>Updating step</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">updating</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-150'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-150'>#</a>
</div>
<p>Retrieve the loglikelihood</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">loglikelihood</span><span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">loglikelihood</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">determinant</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-151'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-151'>#</a>
</div>
<p>Prediction step</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">prediction</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-152'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-152'>#</a>
</div>
<p>Aids to numerical stability</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">numerical_stability</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-153'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-153'>#</a>
</div>
<p>Check for convergence</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">check_convergence</span><span class="p">()</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-154'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-154'>#</a>
</div>
<p>Advance the time</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">t</span> <span class="o">+=</span> <span class="mf">1</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">initialize_object_pointers</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span> <span class="k">except</span> <span class="o">*</span><span class="p">:</span>
<span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">t</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-155'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-155'>#</a>
</div>
<p>Indices for possibly time-varying arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">design_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">obs_intercept_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">obs_cov_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">transition_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">state_intercept_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">selection_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">state_cov_t</span> <span class="o">=</span> <span class="mf">0</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-156'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-156'>#</a>
</div>
<p>Get indices for possibly time-varying arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">time_invariant</span><span class="p">:</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">design</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">design_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">obs_intercept_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">obs_cov_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">transition</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">transition_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">state_intercept</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">state_intercept_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">selection</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">selection_t</span> <span class="o">=</span> <span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">state_cov_t</span> <span class="o">=</span> <span class="n">t</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-157'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-157'>#</a>
</div>
<p>Initialize object-level pointers to statespace arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">obs</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">design</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">design</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">design_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_intercept</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_intercept</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">obs_intercept_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">obs_cov_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">transition</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">transition</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">transition_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_intercept</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">state_intercept</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">state_intercept_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selection</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">selection</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">selection_t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">state_cov_t</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-158'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-158'>#</a>
</div>
<p>Initialize object-level pointers to initialization</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">initialized</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">(</span><span class="s">&quot;Statespace model not initialized.&quot;</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">initial_state</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">initial_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">initial_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-159'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-159'>#</a>
</div>
<p>Initialize object-level pointers to input arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">input_state</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">input_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-160'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-160'>#</a>
</div>
<p>Copy initialization arrays to input arrays if we're starting the
filter</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="n">t</span> <span class="o">==</span> <span class="mf">0</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-161'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-161'>#</a>
</div>
<p><code>predicted_state[:,0]</code> $= a_1 =$ <code>initial_state</code><br />
<code>predicted_state_cov[:,:,0]</code> $= P_1 =$ <code>initial_state_cov</code> </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">initial_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">input_state</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">initial_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">input_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-162'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-162'>#</a>
</div>
<p>Initialize object-level pointers to output arrays</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">forecast</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">forecast</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">forecast_error</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">filtered_state</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">filtered_state</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">filtered_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">predicted_state</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">predicted_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">initialize_function_pointers</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span> <span class="k">except</span> <span class="o">*</span><span class="p">:</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">filter_method</span> <span class="o">&amp;</span> <span class="n">FILTER_CONVENTIONAL</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">forecasting</span> <span class="o">=</span> <span class="n">dforecast_conventional</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">&amp;</span> <span class="n">INVERT_UNIVARIATE</span> <span class="ow">and</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">==</span> <span class="mf">1</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion</span> <span class="o">=</span> <span class="n">dinverse_univariate</span>
<span class="k">elif</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">&amp;</span> <span class="n">SOLVE_CHOLESKY</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion</span> <span class="o">=</span> <span class="n">dsolve_cholesky</span>
<span class="k">elif</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">&amp;</span> <span class="n">SOLVE_LU</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion</span> <span class="o">=</span> <span class="n">dsolve_lu</span>
<span class="k">elif</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">&amp;</span> <span class="n">INVERT_CHOLESKY</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion</span> <span class="o">=</span> <span class="n">dinverse_cholesky</span>
<span class="k">elif</span> <span class="bp">self</span><span class="o">.</span><span class="n">inversion_method</span> <span class="o">&amp;</span> <span class="n">INVERT_LU</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">inversion</span> <span class="o">=</span> <span class="n">dinverse_lu</span>
<span class="k">else</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span><span class="s">&quot;Invalid inversion method&quot;</span><span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">updating</span> <span class="o">=</span> <span class="n">dupdating_conventional</span>
<span class="bp">self</span><span class="o">.</span><span class="n">loglikelihood</span> <span class="o">=</span> <span class="n">dloglikelihood_conventional</span>
<span class="bp">self</span><span class="o">.</span><span class="n">prediction</span> <span class="o">=</span> <span class="n">dprediction_conventional</span>
<span class="k">else</span><span class="p">:</span>
<span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span><span class="s">&quot;Invalid filtering method&quot;</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">select_arrays</span><span class="p">(</span><span class="bp">self</span><span class="p">)</span> <span class="k">except</span> <span class="o">*</span><span class="p">:</span>
<span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">l</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="nb">int</span> <span class="n">design_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">obs_cov_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="nb">int</span> <span class="n">selected_state_cov_t</span> <span class="o">=</span> <span class="mf">0</span>
<span class="n">cdef</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">cdef</span><span class="p">:</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim1</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span>
<span class="n">np</span><span class="o">.</span><span class="n">npy_intp</span> <span class="n">dim2</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-163'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-163'>#</a>
</div>
<h3>Get selected state covariance matrix</h3>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span> <span class="o">==</span> <span class="mf">0</span> <span class="ow">or</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">selected_state_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span>
<span class="n">selected_state_cov_t</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">selected_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">selected_state_cov_t</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-164'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-164'>#</a>
</div>
<h4>Calculate selected state covariance matrix</h4>
<p>$Q_t^* = R_t Q_t R_t'$</p>
<p>Combine the selection matrix and the state covariance matrix to get
the simplified (but possibly singular) "selected" state covariance
matrix (see e.g. Durbin and Koopman p. 43)</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-165'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-165'>#</a>
</div>
<p><code>tmp0</code> array used here, dimension $(m \times r)$ </p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-166'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-166'>#</a>
</div>
<p>$\#_0 = 1.0 * R_t Q_t$<br />
$(m \times r) = (m \times r) (r \times r)$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">selection</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="bp">self</span><span class="o">.</span><span class="n">state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-167'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-167'>#</a>
</div>
<p>$Q_t^* = 1.0 * \#_0 R_t'$<br />
$(m \times m) = (m \times r) (m \times r)'$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dgemm</span><span class="p">(</span><span class="s">&quot;N&quot;</span><span class="p">,</span> <span class="s">&quot;T&quot;</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nposdef</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">alpha</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selection</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="n">beta</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">selected_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">)</span>
<span class="k">else</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_state_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">selected_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-168'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-168'>#</a>
</div>
<h3>Perform missing selections</h3>
<p>In Durbin and Koopman (2012), these are represented as matrix
multiplications, i.e. $Z_t^* = W_t Z_t$ where $W_t$ is a row
selection matrix (it contains a subset of rows of the identity
matrix).</p>
<p>It's more efficient, though, to just copy over the data directly,
which is what is done here. Note that the <code>selected_*</code> arrays are
defined as single-dimensional, so the assignment indexes below are
set such that the arrays can be interpreted by the BLAS and LAPACK
functions as two-dimensional, column-major arrays.</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nmissing</span><span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">0</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-169'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-169'>#</a>
</div>
<p>Set dimensions</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nmissing</span><span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog2</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="o">**</span><span class="mf">2</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendogstates</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">design</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">design_t</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_cov</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mf">2</span><span class="p">]</span> <span class="o">&gt;</span> <span class="mf">1</span><span class="p">:</span> <span class="n">obs_cov_t</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span>
<span class="n">k</span> <span class="o">=</span> <span class="mf">0</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">):</span>
<span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">missing</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_obs</span><span class="p">[</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-170'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-170'>#</a>
</div>
<p>i is rows
k is rows</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nstates</span><span class="p">,</span>
<span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">design</span><span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="n">design_t</span><span class="p">],</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">,</span>
<span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">selected_design</span><span class="p">[</span><span class="n">k</span><span class="p">],</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-171'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-171'>#</a>
</div>
<p>i, k is columns
j, l is rows</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">l</span> <span class="o">=</span> <span class="mf">0</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span><span class="p">):</span>
<span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">missing</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="p">]:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">selected_obs_cov</span><span class="p">[</span><span class="n">l</span> <span class="o">+</span> <span class="n">k</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="p">]</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">obs_cov</span><span class="p">[</span><span class="n">j</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">obs_cov_t</span><span class="p">]</span>
<span class="n">l</span> <span class="o">+=</span> <span class="mf">1</span>
<span class="n">k</span> <span class="o">+=</span> <span class="mf">1</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">selected_obs</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">design</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">selected_design</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="bp">self</span><span class="o">.</span><span class="n">obs_cov</span> <span class="o">=</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">selected_obs_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span>
<span class="k">else</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-172'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-172'>#</a>
</div>
<p>Reset dimensions</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">nendog</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendog2</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span><span class="o">**</span><span class="mf">2</span>
<span class="bp">self</span><span class="o">.</span><span class="n">nendogstates</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">nendog</span> <span class="o">*</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">post_convergence</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-173'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-173'>#</a>
</div>
<p>TODO this should probably be defined separately for each Kalman filter type - e.g. <code>post_convergence_conventional</code>, etc.</p>
</div>
<div class='code'>
<div class="highlight"><pre></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-174'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-174'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-175'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-175'>#</a>
</div>
<p>$F_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nendog2</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">-</span><span class="mf">1</span><span class="p">],</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">forecast_error_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-176'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-176'>#</a>
</div>
<p>$P_{t|t}$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">-</span><span class="mf">1</span><span class="p">],</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">filtered_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-177'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-177'>#</a>
</div>
<p>$P_t$</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span> <span class="mf">0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">-</span><span class="mf">1</span><span class="p">],</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">numerical_stability</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
<span class="k">cdef</span> <span class="kt">int</span> <span class="nf">i</span><span class="p">,</span> <span class="nf">j</span>
<span class="k">cdef</span> <span class="kt">np</span>.<span class="kt">float64_t</span> <span class="nf">value</span>
<span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">stability_method</span> <span class="o">&amp;</span> <span class="n">STABILITY_FORCE_SYMMETRY</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-178'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-178'>#</a>
</div>
<p>Enforce symmetry of predicted covariance matrix<br />
$P_{t+1} = 0.5 * (P_{t+1} + P_{t+1}')$<br />
See Grewal (2001), Section 6.3.1.1</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">):</span>
<span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nstates</span><span class="p">):</span>
<span class="n">value</span> <span class="o">=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="p">(</span>
<span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span> <span class="o">+</span>
<span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span>
<span class="p">)</span>
<span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">value</span>
<span class="bp">self</span><span class="o">.</span><span class="n">model</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="bp">self</span><span class="o">.</span><span class="n">t</span><span class="o">+</span><span class="mf">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">value</span>
<span class="k">cdef</span> <span class="kt">void</span> <span class="nf">check_convergence</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-179'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-179'>#</a>
</div>
<p>Constants</p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">cdef</span><span class="p">:</span>
<span class="nb">int</span> <span class="n">inc</span> <span class="o">=</span> <span class="mf">1</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span>
<span class="n">np</span><span class="o">.</span><span class="n">float64_t</span> <span class="n">gamma</span> <span class="o">=</span> <span class="o">-</span><span class="mf">1.0</span>
<span class="k">if</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">time_invariant</span> <span class="ow">and</span> <span class="ow">not</span> <span class="bp">self</span><span class="o">.</span><span class="n">converged</span><span class="p">:</span></pre></div>
</div>
</div>
<div class='clearall'></div>
<div class='section' id='section-180'>
<div class='docs'>
<div class='octowrap'>
<a class='octothorpe' href='#section-180'>#</a>
</div>
<h4>Check for steady-state convergence</h4>
<p><code>tmp0</code> array used here, dimension $(m \times m)$<br />
<code>tmp1</code> array used here, dimension $(1 \times 1)$ </p>
</div>
<div class='code'>
<div class="highlight"><pre> <span class="n">dcopy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">input_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="n">daxpy</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">gamma</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">predicted_state_cov</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span>
<span class="k">if</span> <span class="n">ddot</span><span class="p">(</span><span class="o">&amp;</span><span class="bp">self</span><span class="o">.</span><span class="n">nstates2</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">tmp0</span><span class="p">,</span> <span class="o">&amp;</span><span class="n">inc</span><span class="p">)</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">tolerance</span><span class="p">:</span>
<span class="bp">self</span><span class="o">.</span><span class="n">converged</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">t</span>
</pre></div>
</div>
</div>
<div class='clearall'></div>
</div>
</body>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment