2a:T473,This notebook illustrates the Kalman filter on one of the simplest state-space models, the **local level model** where the time series has a mean that evolves over time like a random walk: > ${tex` \begin{aligned} y_t &= \mu_t + \varepsilon_t, \hspace{0.75cm} \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \\ \mu_t &= \mu_{t-1} + \eta_t, \quad \eta_t \sim N(0, \sigma_\eta^2) \end{aligned} `} We will simulate data from the local level model for given measurement standard deviation ${tex`\sigma_\varepsilon`} and innovation standard deviation ${tex`\sigma_\eta`}. This gives us a chance to investigate how well we can infer (filter) the true underlying time series of local levels ${tex`\mu_1,\ldots,\mu_T`} from just observing the time series ${tex`y_1,\ldots,y_T`}. In practice we would not know the true parameters ${tex`\sigma_\varepsilon`} and ${tex`\sigma_\eta`} and the widget allows us to see the effect of these choices on the inferred (filtered) local levels ${tex`\mu_1,\ldots,\mu_T`}. Tick the '${tex`\text{estimate }\sigma_\varepsilon \text{ and } \sigma_\eta`}' box to estimate these parameters by maximum likelihood.2b:T4ba,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ₜ", "local level, μₜ", "filter mean", "95 % intervals"], range: [colors[0], colors[1], colors[2], "#F0F0F0", "#d8d8d8"] }, marks: [ Plot.ruleX([0]), 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, tip: false}), Plot.lineY(plotdata, { filter: d => d.type == "local level, μₜ" && plotselector.includes("local level, μₜ"), x: "time", y: "y", stroke: "type", strokeWidth: 2, tip: false}), Plot.lineY(plotdata, { filter: d => d.type == "filter mean" && plotselector.includes("filter mean"), x: "time", y: "y", stroke: "type", strokeWidth: 2, tip: false}), ] })2c:[{"id":"279beaf926475396503df0848853910235ccbcaf8555c1f36196d3a6804a2beba775e04f62a07c6185b8b2ea9300b8e7ac34682a638dfd6b1b2ce505aa189a3c","href":"https://static.observableusercontent.com/files/279beaf926475396503df0848853910235ccbcaf8555c1f36196d3a6804a2beba775e04f62a07c6185b8b2ea9300b8e7ac34682a638dfd6b1b2ce505aa189a3c","name":"nile.csv","type":"text/csv","lastModified":1702895114497.104,"size":929}] 28:{"title":"Kalman filtering simulated local level data","cells":[{"id":0,"value":"# Kalman filtering simulated local level data","pinned":false,"hidden":false,"mode":"md"},{"id":194,"value":"$2a","pinned":false,"hidden":false,"mode":"md"},{"id":463,"value":"#### Settings for the data generating model","pinned":false,"hidden":false,"mode":"md"},{"id":162,"value":"viewof dgp = Inputs.form({\n n: Inputs.range([10, 1000], {value: 100, step: 1, label: \"time series length\"}),\n sigmaeps: Inputs.range([0.01, 5], {value: 1, step: 0.1, label: tex`\\text{true }\\sigma_\\varepsilon`}),\n sigmaeta: Inputs.range([0.01, 2], {value: 1, step: 0.1, label: tex`\\text{true }\\sigma_\\eta`}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1508,"value":"#### Settings for the model used for filtering","pinned":false,"hidden":false,"mode":"md"},{"id":1504,"value":"viewof model = {\n if (optparam){\n return Inputs.form([Inputs.text(\n {value: \"Parameters estimated by maximum likelihood\", label: \"\", disabled: true}\n )])\n }else{\n return Inputs.form({\n sigmaeps: Inputs.range([0.01, 5], \n {value: 1, step: 0.1, label: tex`\\text{fit based on }\\sigma_\\varepsilon`}),\n sigmaeta: Inputs.range([0.01, 5], \n {value: 1, step: 0.1, label: tex`\\text{fit based on }\\sigma_\\eta`}),\n initmean: Inputs.range([-10,10], \n {value: 0, step: 1, label: tex`\\text{initial mean }\\mu_{0|0}`}),\n initstd: Inputs.range([1,1000], \n {value: 100, step: 1, label: tex`\\text{initial stdev }\\sigma_{0|0}`})\n })\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1738,"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ₜ\", \"local level, μₜ\", \"filter mean\", \"95 % intervals\"], {label: \"plot:\", value: [\"time series, xₜ\", \"local level, μₜ\", \"filter mean\", \"95 % intervals\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":781,"value":"viewof simulatebutton = Inputs.button(\"Simulate new data\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1378,"value":"$2b","pinned":false,"hidden":false,"mode":"ojs"},{"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":1774,"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":1728,"value":"import {KalmanFilterUpdate, kalmanFilter} from \"@mattiasvillani/statistics-tools\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":1733,"value":"import {simLocalLevel} from \"@mattiasvillani/statistics-tools\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":96,"value":"#### Data","pinned":false,"hidden":false,"mode":"md"},{"id":1754,"value":"sigmaeps = {\n if (optparam){\n return Math.sqrt(Math.exp(postMode[0]))\n }else{\n return model.sigmaeps\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1756,"value":"sigmaeta = {\n if (optparam){\n return Math.sqrt(Math.exp(postMode[1]))\n }else{\n return model.sigmaeta\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1766,"value":"initmean = {\n if (optparam){\n return 1000\n }else{\n return model.initmean\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1768,"value":"initstd = {\n if (optparam){\n return 1000\n }else{\n return model.initstd\n }\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1760,"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":1708,"value":"colors = [\"#6C8EBF\", \"#c0a34d\", \"#780000\", \"#007878\", \"#bf9d6c\", \"#3A6B35\", \"#b5c6df\",\"#eadaaa\", \"#bb989a\", \"#98bbb9\", \"#bf8d6c\", \"#CBD18F\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1013,"value":"n = dgp.n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1453,"value":"y = simulatedData.filter(d => d.type == \"time series, xₜ\").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: t+1, 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 truelevel = simulatedData.filter(d => d.type == \"local level, μₜ\").map(d => d.x);\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: t+1, y: y[t], type: \"time series, xₜ\"}) \n data.push({time: t+1, y: truelevel[t], type: \"local level, μₜ\"})\n data.push({time: t+1, y: μ_filter[t], type: \"filter mean\"})\n data.push({time: t+1, y: μ_filter[t] - 1.96*std_filter[t], type: \"filter 95% lower\"})\n data.push({time: t+1, 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":519,"value":"simulatedData = {\n simulatebutton\n return simLocalLevel(n, dgp.sigmaeps, dgp.sigmaeta)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1724,"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"},{"id":1762,"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"}],"files":"$W2c","stdlib":"1","titling":"auto"}