2a:T420,plotfilter = Plot.plot({ title: optparam ? md`Maximum likelihood estimates: ${tex`\hat\sigma_\varepsilon = ${sigmaeps.toFixed(2)}`} and ${tex`\hat\sigma_\eta = ${sigmaeta.toFixed(2)}`}` : "", color: { legend: true, domain: ["time series, xₜ", "filter mean", "95 % intervals"], range: [colors[0], colors[2], "#d8d8d8"] }, y: { label: "river flow in 10⁸m³" }, x: { label: "year" }, marks: [ Plot.ruleX([new Date(1870 , 0, 2)]), Plot.ruleY([d3.min(plotdata.map(d => d.y))]), Plot.areaY(plotdataInterval, {filter: plotselector.includes("95 % intervals"), x: "time", y1: "lower", y2: "upper", fill: "#d8d8d8"}), Plot.lineY(plotdata, { filter: d => d.type == "time series, xₜ" && plotselector.includes("time series, xₜ"), x: "time", y: "y", stroke: "type", strokeWidth: 2}), Plot.lineY(plotdata, { filter: d => d.type == "filter mean" && plotselector.includes("filter mean"), x: "time", y: "y", stroke: "type", strokeWidth: 2, tip: false}), ] })2b:[{"id":"279beaf926475396503df0848853910235ccbcaf8555c1f36196d3a6804a2beba775e04f62a07c6185b8b2ea9300b8e7ac34682a638dfd6b1b2ce505aa189a3c","href":"https://static.observableusercontent.com/files/279beaf926475396503df0848853910235ccbcaf8555c1f36196d3a6804a2beba775e04f62a07c6185b8b2ea9300b8e7ac34682a638dfd6b1b2ce505aa189a3c","name":"nile.csv","type":"text/csv","lastModified":1702897538444.16,"size":929}] 28:{"title":"Kalman filtering Nile river data","cells":[{"id":0,"value":"# Kalman filtering Nile river data","pinned":false,"hidden":false,"mode":"md"},{"id":1789,"value":"### Interactive exploration of the Kalman filter","pinned":false,"hidden":false,"mode":"md"},{"id":194,"value":"We will here illustrate the Kalman filter on the Nile river data using the **local level model** where the time series has a mean that evolves over time like a random walk:\n\n> ${tex`\n\\begin{aligned}\n y_t &= \\mu_t + \\varepsilon_t, \\hspace{0.75cm} \\varepsilon_t \\sim N(0, \\sigma_\\varepsilon^2) \\\\\n \\mu_t &= \\mu_{t-1} + \\eta_t, \\quad \\eta_t \\sim N(0, \\sigma_\\eta^2)\n\\end{aligned}\n`}\n \nThe Kalman filter (see below for more details) uses the observed data ${tex`y_1,\\ldots,y_t`} to estimate the local mean ${tex`\\mu_t`} at time *t*. Note that this is an exercise in **filtering**, which is the *instantaneous* estimate at time *t*, using only the data available up to that time period. ","pinned":false,"hidden":false,"mode":"md"},{"id":1777,"value":"Here are few things to try out in the widget below:\n- Increase the measurement standard deviation ${tex`\\sigma_\\varepsilon`} and note how the filtered mean pays less and less attention to changes in the data (because the model believes that the data is very poor quality (noisy) and tells us basically nothing about the level). Then move ${tex`\\sigma_\\varepsilon`} to smaller values and note how the filtered mean starts chasing the data (because the model believes that the data are super informative about the level).\n- Make the standard deviation for the initial level ${tex`\\sigma_0`} very small and then change the initial mean ${tex`\\mu_0`} to see how this affects the filtered mean at the first part of the time series.\n- Move the standard deviation of the innovations to the level ${tex`\\sigma_\\eta`} small and note how the filtered mean becomes smoother and smoother over time.\n- Tick the '${tex`\\text{estimate }\\sigma_\\varepsilon \\text{ and } \\sigma_\\eta`}' box to estimate these parameters by maximum likelihood.","pinned":false,"hidden":false,"mode":"md"},{"id":1508,"value":"##### Model settings","pinned":false,"hidden":false,"mode":"md"},{"id":1983,"value":"viewof param_input = {\n if (optparam){\n return Inputs.form([Inputs.text({value: \"Parameters estimated by maximum likelihood\", label: \"\", disabled: true})])\n }else{\n return Inputs.form({\n sigmaeps: Inputs.range([1, 300], {value: 100, step: 0.1, label: tex`\\text{noise std }\\sigma_\\varepsilon`}),\n sigmaeta: Inputs.range([1, 300], {value: 100, step: 0.1, label: tex`\\text{innov std }\\sigma_\\eta`}),\n initmean: Inputs.range([-2000,2000], {value: 1000, step: 1, label: tex`\\text{initial mean }\\mu_{0|0}`}),\n initstd: Inputs.range([1,1000], {value: 100, step: 1, label: tex`\\text{initial stdev }\\sigma_{0|0}`})\n })\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1957,"value":"viewof optparam = Inputs.toggle({label: tex`\\text{estimate }\\sigma_\\varepsilon \\text{ and } \\sigma_\\eta`, value: false})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1531,"value":"viewof plotselector = Inputs.checkbox([\"time series, xₜ\", \"filter mean\", \"95 % intervals\"], {label: \"plot:\", value: [\"time series, xₜ\", \"filter mean\", \"95 % intervals\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1378,"value":"$2a","pinned":false,"hidden":false,"mode":"ojs"},{"id":1794,"value":"#### State-space models and the Kalman filter","pinned":false,"hidden":false,"mode":"md"},{"id":1189,"value":"There are two different types of relevant inferences for the unobserved state ${tex`\\mu_t`} in the local level model:\n\n- **Filtering** which is the **instantaneous** estimate of ${tex`\\mu_t`} given all the data ${tex`y_1,\\ldots,y_t`} up to time *t*. We denote this estimate by ${tex`\\mu_{t|t}`} to show that\n - it is the inference for ${tex`\\mu`} at time *t* (the index before the | sign),\n - using the data up to time *t* (the index after the | sign).\n- **Smoothing** which is the **retrospective** estimate of ${tex`\\mu_t`} given all the data ${tex`y_1,\\ldots,y_T`} until the end of the time series *T*. We denote this estimate by ${tex`\\mu_{t|T}`}.","pinned":false,"hidden":false,"mode":"md"},{"id":1200,"value":"The local level model is an example of a Linear Gaussian State-Space (LGSS) model\n\n> ${tex`\n\\begin{aligned}\ny_t &= \\boldsymbol{F} \\boldsymbol{\\theta}_t + v_t,\\hspace{1.3cm} v_t \\sim N(\\boldsymbol{0},\\boldsymbol{V}) \\\\\n\\boldsymbol{\\theta}_t &= \\boldsymbol{G} \\boldsymbol{\\theta}_{t-1} + \\boldsymbol{w}_t, \\qquad \\boldsymbol{w}_t \\sim N(\\boldsymbol{0},\\boldsymbol{W})\n\\end{aligned}\n`}\n\nFor example, the local level model is a state-space model with\n\n> ${tex`\n\\boldsymbol{\\theta}_t = \\mu_t, \n \\boldsymbol{F} = \\boldsymbol{G} = 1,\n \\boldsymbol{V} = \\sigma_\\varepsilon^2,\n\\text{ and }\n\\boldsymbol{W} = \\sigma_\\nu^2.\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":1207,"value":"There are three main inferences of interest in the state-space model:\n- The **filtered** estimate ${tex`\\hat{\\boldsymbol{\\theta}}_{t|t}`} of the state ${tex`\\boldsymbol{\\theta}_t`} uses **data up to time** *t*.\n- The **smoothed** estimate ${tex`\\hat{\\boldsymbol{\\theta}}_{t|T}`} of the state ${tex`\\boldsymbol{\\theta}_t`} uses **all the data** until the end of the time series *T*.\n- Inference for the parameters in the state-space model, i.e. the parts of ${tex`\\boldsymbol{F}, \\boldsymbol{G},\n \\boldsymbol{V} \\text{ and } \\boldsymbol{W}`} which are unknown and need to be estimated.","pinned":false,"hidden":false,"mode":"md"},{"id":1215,"value":"We will here be concerned with **filtering**, i.e. computing the sequence of instantaneous estimates of the unobserved state at every time point ${tex`t=1,2,\\ldots,T`}:\n\n> ${tex`\\hat{\\boldsymbol{\\theta}}_{1|1},\\hat{\\boldsymbol{\\theta}}_{2|2},\\ldots,\\hat{\\boldsymbol{\\theta}}_{T|T}.`}\n\nFor the linear Gaussian state-space model, the filtering can be efficiently computed by the **Kalman filter**. We will use the notation ${tex`\\boldsymbol{\\mu}_{t|t}`} to denote the point estimates of the state vector (posterior mean in a Bayesian setting) and ${tex`\\boldsymbol{\\Omega}_{t|t}`} for the covariance of the state vector (posterior covariance in a Bayesian setting).
\nThe Kalman filter algorithm computes the estimates recursively in time, starting with some initial guess ${tex`\\boldsymbol{\\mu}_{0|0}`} of the state and some uncertainty/covariance matrix ${tex`\\boldsymbol{\\Omega}_{0|0}`} for that guess. Here is the Kalman filter algorithm in full:","pinned":false,"hidden":false,"mode":"md"},{"id":1253,"value":"---\n**Initialization**: set ${tex`\\boldsymbol{\\mu}_{0|0}`} and ${tex`\\boldsymbol{\\Omega}_{0|0}`} \n\n`for` ${tex`t=1,\\ldots,T`} `do`\n\n> **Prediction update**\n\n> ${tex`\n\\begin{aligned}\n\\boldsymbol{\\mu}_{t|t-1} &= \\boldsymbol{G} \\boldsymbol{\\mu}_{t-1|t-1} \\\\ \n\\boldsymbol{\\Omega}_{t|t-1} &= \\boldsymbol{G}\\boldsymbol{\\Omega}_{t-1|t-1} \\boldsymbol{G}^\\top + \\boldsymbol{W}\n\\end{aligned}\n`}\n\n> **Measurement update**\n\n> ${tex`\n\\begin{aligned}\n\\boldsymbol{\\mu}_{t|t} &= \\boldsymbol{\\mu}_{t|t-1} + \\boldsymbol{K}_t ( y_t - \\boldsymbol{F} \\boldsymbol{\\mu}_{t|t-1} ) \\\\ \n\\boldsymbol{\\Omega}_{t|t} &= (\\boldsymbol{I} - \\boldsymbol{K}_t \\boldsymbol{F} )\\boldsymbol{\\Omega}_{t|t-1} \n\\end{aligned}\n`}\n\n> where ${tex`\\boldsymbol{K}_t = \\boldsymbol{\\Omega}_{t|t-1}\\boldsymbol{F}^\\top ( \\boldsymbol{F} \\boldsymbol{\\Omega}_{t|t-1}\\boldsymbol{F}^\\top + \\boldsymbol{V})^{-1}`} is the **Kalman Gain**.\n\n`end`\n\n---","pinned":false,"hidden":false,"mode":"md"},{"id":13,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\nHere follow code and details that you probably do not want to read.","pinned":false,"hidden":false,"mode":"md"},{"id":15,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":11,"value":"math = require(\"mathjs\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":816,"value":"jStat = require(\"jstat\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":1861,"value":"optimjs = require('optimization-js@1.5.0/dist/optimization.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1231,"value":"import {textcolor} from \"@observablehq/text-color-annotations-in-markdown\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2120,"value":"import {KalmanFilterUpdate, kalmanFilter} from \"@mattiasvillani/statistics-tools\"\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":96,"value":"#### Data","pinned":false,"hidden":false,"mode":"md"},{"id":1735,"value":"colors = [\"#6C8EBF\", \"#c0a34d\", \"#780000\", \"#007878\", \"#bf9d6c\", \"#3A6B35\", \"#b5c6df\",\"#eadaaa\", \"#bb989a\", \"#98bbb9\", \"#bf8d6c\", \"#CBD18F\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1013,"value":"n = y.length","pinned":false,"hidden":false,"mode":"ojs"},{"id":1711,"value":"data = FileAttachment(\"nile.csv\").csv({typed: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1709,"value":"y = data.map(d => d.x)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1311,"value":"filter_results = {\nconst Y = y; // Your measurement data\nconst A = math.identity(1); // Your system dynamics matrix\nconst C = math.identity(1); // Your measurement matrix\nconst Σₑ = math.matrix([sigmaeps**2]) // Process noise covariance matrix\nconst Σₙ = math.matrix([sigmaeta**2]); // Measurement noise covariance matrix\nconst μ0 = math.matrix([initmean]); // Initial mean state vector\nconst Σ0 = math.matrix([initstd**2]); // Initial state covariance matrix\n\nconst { μ_filter, Σ_filter } = kalmanFilter(Y, A, C, Σₑ, Σₙ, μ0, Σ0);\n return [μ_filter, Σ_filter]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1484,"value":"plotdataInterval = {\n const data = [];\n const lower = plotdata.filter(d => d.type == \"filter 95% lower\").map(d => d.y);\n const upper = plotdata.filter(d => d.type == \"filter 95% upper\").map(d => d.y);\n for (let t = 0; t < n; ++t) {\n data.push({time: new Date(1871 + t, 0, 2), lower: lower[t], upper: upper[t]})\n }\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1348,"value":"plotdata = {\n var data = [];\n const μ_filter = filter_results[0].map(d => d._data[0]);\n const std_filter = filter_results[1].map(d => Math.sqrt(d._data[0]));\n for (let t = 0; t < n; ++t) {\n data.push({time: new Date(1871 + t, 0, 2) , y: y[t], type: \"time series, xₜ\"}) \n data.push({time: new Date(1871 + t, 0, 2), y: μ_filter[t], type: \"filter mean\"})\n data.push({time: new Date(1871 + t, 0, 2), y: μ_filter[t] - 1.96*std_filter[t], type: \"filter 95% lower\"})\n data.push({time: new Date(1871 + t, 0, 2), y: μ_filter[t] + 1.96*std_filter[t], type: \"filter 95% upper\"})\n }\n return data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":94,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md"},{"id":1895,"value":"function negLogLik(param){\n \n const Y = y; // Your measurement data\n const A = math.identity(1); // Your system dynamics matrix\n const C = math.identity(1); // Your measurement matrix\n const Σₑ = math.matrix([Math.exp(param[0])]) // Measurement noise covariance matrix\n const Σₙ = math.matrix([Math.exp(param[1])]); // Process noise covariance matrix \n const μ0 = math.matrix([initmean]); // Initial mean state vector\n const Σ0 = math.matrix([initstd**2]); // Initial state covariance matrix\n \n const { μ_filter, Σ_filter, μ_prop, Σ_prop } = kalmanFilter(Y, A, C, Σₑ, Σₙ, μ0, Σ0);\n\n const CΣₑC = math.multiply( math.multiply(C, Σₑ), math.transpose(C) );\n let logLik = 0;\n for (let t = 0; t < y.length; ++t) {\n logLik += normal_logpdf( y[t], μ_prop[t]._data[0], Math.sqrt( math.add( Σ_prop[t], CΣₑC )._data[0] ) )\n }\n return -logLik\n \n }","pinned":false,"hidden":false,"mode":"ojs"},{"id":1901,"value":"postMode = {\n if (optparam){\n return optimjs.minimize_Powell(negLogLik, [10,10]).argument\n }else{\n return NaN\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2035,"value":"initmean = {\n if (optparam){\n return 1000\n }else{\n return param_input.initmean\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2039,"value":"initstd = {\n if (optparam){\n return 1000\n }else{\n return param_input.initstd\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1917,"value":"sigmaeps = {\n if (optparam){\n return Math.sqrt(Math.exp(postMode[0]))\n }else{\n return param_input.sigmaeps\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1971,"value":"sigmaeta = {\n if (optparam){\n return Math.sqrt(Math.exp(postMode[1]))\n }else{\n return param_input.sigmaeta\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1886,"value":"function normal_logpdf(x, mu, sigma){\n return -0.5*Math.log(2*Math.PI) - Math.log(sigma) -0.5*( ((x-mu)/sigma)**2 )\n}","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W2b","stdlib":"1","titling":"auto"}