2b:null 2c:["$","title",null,{"children":"Autoregressive processes - simulation and priors | Mattias Villani | Observable"}] 31:T421,viewof phi = { switch (siminputs[0]) { case 1: return Inputs.form([ Inputs.range([-2, 2], {value: 0.7, step: 0.01, label: tex`\phi_1`}) ]) break; case 2: return Inputs.form([ Inputs.range([-2, 2], {value: 0.7, step: 0.01, label: tex`\phi_1`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_2`}) ]) break; case 3: return Inputs.form([ Inputs.range([-2, 2], {value: 0.7, step: 0.01, label: tex`\phi_1`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_2`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_3`}) ]) break; case 4: return Inputs.form([ Inputs.range([-2, 2], {value: 0.7, step: 0.01, label: tex`\phi_1`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_2`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_3`}), Inputs.range([-2, 2], {value: 0.0, step: 0.01, label: tex`\phi_4`}) ]) break; } } 32:T448,## Prior We need to place a prior on all model parameters: ${tex`\mu`}, ${tex`\phi_1,\ldots,\phi_p`} and ${tex`\sigma^2`}.
A commonly used prior assumes independence between all parameters and: - a subjective prior on the unconditional mean ${tex.block`\mu \sim N(m_0,s_0^2)`} from expert information. If the process is stationary, then the long-range forecasts will converge to ${tex`\mu`}, which can be a useful way to think about ${tex`\mu`} for the expert. - a prior for the standard deviation of the errors ${tex`\varepsilon`}, for example the scaled inverse chi-squared distribution ${tex.block`\sigma^2 \sim \mathrm{Inv-}\chi^2(\nu_0,\sigma_0^2)`} often with a low degrees of freedom ${tex`\nu_0`} (so that the prior is rather vague) - a regularization prior on the AR coefficients ${tex`\phi_1,\ldots,\phi_p`} that has the effect of shrinking more towards zero for the longer lags, i.e. on ${tex`\phi_k`} for larger *k*: ${tex.block`\begin{aligned} \phi_1 &\sim N(\mu_0,\sigma^2\lambda^2)\\ \phi_k &\sim N\Big(0,\sigma^2\frac{\lambda^2}{k^\gamma}\Big) \text{ for }k>1 \end{aligned}`}33:[] 2f:{"id":"112b4245ebf148d9","slug":"ar-processes","title":"Autoregressive processes - simulation and priors","thumbnail":"e25fcb72a471b447a1c486e9fe813ad905a50609b6130a9ed75dc87a0718f95d","access_level":"public","create_time":"$D2023-04-06T13:09:43.819Z","update_time":"$D2025-02-09T20:31:38.432Z","view_time":null,"first_public_time":"$D2023-04-06T13:09:43.819Z","fork":{"id":"4594bdcee9954fb1","slug":"student-t-distribution","title":"Student-t distribution (three-parameter)","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"stars":0,"starred":false,"owner":{"type":"team","id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":"pro_2024"},"version":939,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null}],"members":[{"type":"individual","id":"56499be8cf43f0ef","login":"mattiasvillani","name":"Mattias Villani","avatar_url":"https://avatars-next.observableusercontent.com/a440aaf39a3abefb86cbcc758fd7a5ab42555ae68d6825b52ec006f385a94d1d","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":2,"collection":{"id":"9fdb7c9ba7be6ba8","slug":"bayesian-learning","title":"Bayesian Learning","owner":{"id":"7592a7ab121f74a0","login":"mattiasvillani","name":"Mattias Villani"}},"type":"notebook","body":{"title":"Autoregressive processes - simulation and priors","cells":[{"id":38,"value":"# Autoregressive processes - simulation and priors","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":406,"value":"## Model\nLet the model be an autoregressive process of order *p* in so called steady-state form where ${tex`\\mu`} is the unconditional mean of the process:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":408,"value":"model_text = tex`\n\\begin{aligned}\nX_t = \\mu + \\phi_1(X_{t-1}-\\mu)+ \\ldots + \\phi_p(X_{t-p}-\\mu) + \\varepsilon_t, \\, \\text{ where } \\, \\varepsilon_t \\overset{\\mathrm{iid}}{\\sim}N(0,\\sigma^2)\n\\end{aligned}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":668,"value":"Here is a simulator for the AR(p) process:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":723,"value":"viewof T = Inputs.range([5, 1000], {value: 50, step: 1, label: tex`\\text{length series, }T`})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":666,"value":"viewof siminputs = Inputs.form([\n Inputs.range([1, 4], {value: 3, step: 1, label: tex`\\text{number of lags, }p`}),\n Inputs.range([-3, 3], {value: 0.0, step: 0.1, label: tex`\\text{mean, }\\mu`}),\n Inputs.range([0, 2], {value: 0.5, step: 0.1, label: tex`\\text{error stdev, }\\sigma`})\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":672,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":784,"value":"\n \n \n \n

\n ${plt_timeseries}\n

\n \n

\n ${pltacf}\n

\n \n","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined"},{"id":521,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":929,"value":"where ${tex`\\lambda>0`} determines the standard deviation of the coefficient on the first lag, ${tex`\\gamma>0`} is hyperparameter that determines how quickly the variance decreases with the lag. Note how the mean for the first lag is allowed to be non-zero, but higher order lags are all centered around zero. This corresponds to the prior belief that the most likely process is the AR(1):\n${tex.block`X_t = \\mu + \\mu_0 (X_{t-1}-\\mu) + \\varepsilon_t, \\, \\text{ where } \\, \\varepsilon_t \\overset{\\mathrm{iid}}{\\sim}N(0,\\sigma^2)`}","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":618,"value":"You can experiment with the prior hyperparameters below to see the effect on the prior.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":118,"value":"viewof myinputs = Inputs.form([\n Inputs.range([1, 10], {value: 3, step: 1, label: tex`\\text{number of lags, }p`}),\n Inputs.range([-1, 1], {value: 0.5, step: 0.1, label: tex`\\text{mean first lag, }\\mu_0`}),\n Inputs.range([0, 2], {value: 0.5, step: 0.1, label: tex`\\text{stdev first lag, }\\lambda`}),\n Inputs.range([0, 5], {value: 2, step: 0.1, label: tex`\\text{lag decay }\\gamma`})\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":305,"value":"plt = Plot.plot({\n color: {\n legend: true,\n type: \"categorical\",\n //domain: d3.range(10).map((d) => `Category ${d + 1}`),\n //scheme: \"tableau10\" \n },\n x: {\n label: \"phi\",\n axis: true\n },\n y: {\n axis: false,\n //domain: [0,1.05*jstat.normal.pdf(myinputs[0], myinputs[0], myinputs[1])]\n },\n marks: [\n Plot.ruleY([0]),\n Plot.line(pdfdata, {x: \"phi\", y: \"priordens\", stroke : \"lag\", strokeWidth: 2}),\n ]\n })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":882,"value":"## Prior predictive draws\nThe graph below simulates time series by first drawing from the prior and the drawing a time series of a certain length from the generated parameter draw.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":905,"value":"viewof Tprior = Inputs.range([10, 1000], {value: 200, step: 1, label: \"time series length for prior predictive\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":834,"value":"viewof priorpredinputs = Inputs.form([\n Inputs.range([0.01, 2], {value: 0.5, step: 0.1, label: tex`\\text{location error std, }\\sigma_0`}),\n Inputs.range([1, 30], {value: 3, step: 1, label: tex`\\text{df error std, }\\nu_0`}),\n Inputs.range([-5, 5], {value: 1, step: 0.1, label: tex`\\text{mean of mean, }m_0`}),\n Inputs.range([0.01, 5], {value: 1, step: 0.1, label: tex`\\text{std of mean, }s_0`})\n ])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":862,"value":"viewof simnow = Inputs.button(\"Simulate new realization!\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":817,"value":"plt_priorpred_timeseries = Plot.plot({\n title: \"Prior predictive realization\",\n x: {\n label: \"time, t\",\n axis: true, \n domain: [0,Tprior]\n },\n y: {\n label: \"X(t)\",\n axis: true,\n },\n marks: [\n //Plot.ruleY([d3.min(timeseries.map(d => d.x))]),\n Plot.ruleX([0]),\n Plot.lineY(timeseries_priorpred, {filter: d => d.time < Tprior, x: \"time\", y: \"x\", stroke: \"steelblue\", strokeWidth: 2}),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":109,"value":"
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n---\n## 🙈 Appendix 🙈\n\nHere is only code and details that you probably don't want to read.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":792,"value":"#### Plots","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":694,"value":"plt_timeseries = Plot.plot({\n x: {\n label: \"time, t\",\n axis: true, \n domain: [0,T]\n },\n y: {\n label: \"X(t)\",\n axis: true,\n },\n marks: [\n Plot.ruleY([d3.min(timeseries.map(d => d.x))]),\n Plot.ruleX([0]),\n Plot.lineY(timeseries, {filter: d => d.time < T, x: \"time\", y: \"x\", stroke: \"steelblue\", strokeWidth: 2})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":771,"value":"pltacf = Plot.plot({\n //title: md`Estimated ACF of AR(${p}) process`,\n marks: [\n //Plot.axisX({ticks: d3.range(n_acf_lags)}),\n Plot.ruleY([0]),\n Plot.ruleX([0]),\n //Plot.lineY(acf.map((d,i) => ({lag: i, acf: d})), {x: \"lag\", y: \"acf\", stroke: \"green\"}),\n Plot.dot(acf.map((d,i) => ({lag: i, acf: d})), {filter: d => d.lag > 0, x: \"lag\", y: \"acf\", fill: \"green\", stroke: \"green\"}),\n Plot.ruleX(acf.map((d,i) => ({lag: i, acf: d})), {filter: d => d.lag > 0, x: \"lag\", y: \"acf\", stroke: \"green\", strokeWidth: 2}),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":221,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":26,"value":"jstat = require('jstat')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":747,"value":"mathjs = require('mathjs')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":778,"value":"import {computeACF, simulateAR} from \"@mattiasvillani/statistics-tools\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":620,"value":"### Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":223,"value":"#### Data","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":547,"value":"p = myinputs[0]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":549,"value":"mu0 = myinputs[1]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":551,"value":"tau = myinputs[2]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":553,"value":"gamma = myinputs[3]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":689,"value":"mu = siminputs[1]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":684,"value":"sigma = siminputs[2]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":895,"value":"sigma0 = priorpredinputs[0]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":897,"value":"nu0 = priorpredinputs[1]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":899,"value":"m0 = priorpredinputs[2]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":901,"value":"s0 = priorpredinputs[3]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":662,"value":"timeseries = simulateAR(1000, mu, phi, sigma)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":17,"value":"pdfdata = {\n var pdfs;\n const phi_grid = d3.range(-3, 3, 0.01);\n pdfs = phi_grid.map(phi => ({phi: phi, priordens: jstat.normal.pdf(phi, mu0, tau), lag: \"lag \"+1}));\n var tmp;\n for (let k = 2; k <= p; k++){\n tmp = phi_grid.map(phi => ({phi: phi, priordens: jstat.normal.pdf(phi, 0, Math.sqrt(tau**2/(k**gamma))), lag: \"lag \"+ k}))\n pdfs = pdfs.concat(tmp)\n }\n return pdfs\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":775,"value":"n_acf_lags = Math.max(5*p,10)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":773,"value":"acf = computeACF(timeseries.map(d => d.x), n_acf_lags)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":822,"value":"timeseries_priorpred = {\n simnow;\n \n let phi = [];\n phi.push(jstat.normal.sample(mu0, tau))\n for (let k = 2; k<=p; k++){\n phi.push( jstat.normal.sample(0, Math.sqrt(tau**2/(k**gamma))))\n }\n let sigma = Math.sqrt(nu0*(sigma0**2)*(1/jstat.chisquare.sample(nu0)))\n let mu = jstat.normal.sample(m0,s0);\n let x = simulateAR(1000, mu, phi, sigma)\n \n return x\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":229,"value":"#### Style","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined"},{"id":202,"value":"","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined"},{"id":58,"value":"viewof md_colour = Inputs.color({label: \"Highlight color cells\", value: \"#F5F5DC\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 30:"$2f:body"