2a:T487,viewof ϕvect = { switch (input.p) { case 0: return Inputs.form([Inputs.text({value: "No regular AR", label: "", disabled: true})]) case 1: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\phi_1`}) ]) break; case 2: return Inputs.form([ Inputs.range([-2, 2], {value: 0.5, step: 0.01, label: tex`\phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_2`}) ]) break; case 3: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_3`}) ]) break; case 4: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_3`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\phi_4`}) ]) break; } }2b:T49b,viewof θvect = { switch (input.q) { case 0: return Inputs.form([Inputs.text({value: "No regular MA", label: "", disabled: true})]) case 1: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\theta_1`}) ]) break; case 2: return Inputs.form([ Inputs.range([-2, 2], {value: 0.5, step: 0.01, label: tex`\theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_2`}) ]) break; case 3: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_3`}) ]) break; case 4: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_3`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\theta_4`}) ]) break; } }2c:T488,viewof Φvect = { switch (input.P) { case 0: return Inputs.form([Inputs.text({value: "No seasonal AR", label: "", disabled: true})]) case 1: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Phi_1`}) ]) break; case 2: return Inputs.form([ Inputs.range([-2, 2], {value: 0.5, step: 0.01, label: tex`\Phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_2`}) ]) break; case 3: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_3`}) ]) break; case 4: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Phi_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_3`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Phi_4`}) ]) break; } }2d:T49c,viewof Θvect = { switch (input.Q) { case 0: return Inputs.form([Inputs.text({value: "No seasonal MA", label: "", disabled: true})]) case 1: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Theta_1`}) ]) break; case 2: return Inputs.form([ Inputs.range([-2, 2], {value: 0.5, step: 0.01, label: tex`\Theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_2`}) ]) break; case 3: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_3`}) ]) break; case 4: return Inputs.form([ Inputs.range([-1, 1], {value: 0.5, step: 0.01, label: tex`\Theta_1`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_2`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_3`}), Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: tex`\Theta_4`}) ]) break; } }2e:T453,function SpecDens(ω, ϕvect, θvect, Φvect, Θvect, σ, s){ const eiw = math.exp(math.complex(0,-ω)) const eiws = math.exp(math.complex(0,-s*ω)) // Regular AR var ARpoly = math.complex(1,0); var lag = 0; for (let ϕ of ϕvect){ lag = lag + 1; ARpoly = math.add(ARpoly, math.multiply(math.complex(-ϕ,0),math.pow(eiw,lag))) } // Seasonal AR var ARspoly = math.complex(1,0); var lag = 0; for (let Φ of Φvect){ lag = lag + 1; ARspoly = math.add(ARspoly, math.multiply(math.complex(-Φ,0),math.pow(eiws,lag))) } // Regular MA var MApoly = math.complex(1,0); var lag = 0; for (let θ of θvect){ lag = lag + 1; MApoly = math.add(MApoly, math.multiply(math.complex(θ,0),math.pow(eiw,lag))) } // Seasonal MA var MAspoly = math.complex(1,0); var lag = 0; for (let Θ of Θvect){ lag = lag + 1; MAspoly = math.add(MAspoly, math.multiply(math.complex(Θ,0),math.pow(eiws,lag))) } return (σ**2/Math.PI)* ((math.abs(MApoly)**2 * math.abs(MAspoly)**2) / (math.abs(ARpoly)**2 * math.abs(ARspoly)**2) ) }2f:T482,function simARMA(n, ϕvect, θvect, Φvect, Θvect, mu, σ) { simulatebutton // Compute product of regular and seasonal polynomials, including all zero coeffs const arparam = combineARPoly(ϕvect, Φvect).map(d => -d) const maparam = combineMAPoly(θvect, Θvect) var p = arparam.length var q = maparam.length const timeSeries = []; const maxpq = Math.max(p,q) // Generate white noise sequence const ε = d3.range(maxpq + n).map(d => σ*d3.randomNormal()()) // Initialize the first p values of the time series for (const j of d3.range(maxpq)){ timeSeries[j] = {time: j - maxpq + 1, x: mu}; } for (let t = maxpq; t < n + maxpq; t++) { let arComponent = 0 if (p > 0){ for (const [index, ϕ] of arparam.entries()){ arComponent = arComponent + ϕ * timeSeries[t - index - 1].x; } } let maComponent = 0 if (q > 0){ for (const [index, θ] of maparam.entries()){ maComponent = maComponent + θ * ε[t - index - 1]; } } timeSeries[t] = {time: t - maxpq + 1, x: mu + arComponent + maComponent + ε[t]} ; } return timeSeries.slice(p); }30:[] 28:{"title":"Spectral density of a seasonal ARMA(p,q)×(P,Q) process","cells":[{"id":0,"value":"# Spectral density of a seasonal ARMA(p,q)×(P,Q) process","pinned":false,"hidden":false,"mode":"md"},{"id":194,"value":"The seasonal ARMA(p,q)×(P,Q) process with seasonality *s* is of the form\n\n${tex`\n\\phi_p(L)\\Phi_P(L)(y_t-\\mu) = \\theta_q(L)\\Theta_Q(L)\\varepsilon_t,\n`}\n\nwhere
\n${tex`\\phi_p(L) = 1 - \\phi_1L - \\cdots - \\phi_pL^p`} is the AR lag polynomial,
\n${tex`\\theta_q(L) = 1 + \\theta_1L + \\cdots + \\theta_qL^q`}is the MA \nlag polynomial
\n${tex`\\Phi_P(L) = 1 - \\Phi_1 L^s - \\cdots - \\Phi_P L^{sP}`} is the seasonal AR lag polynomial,
\n${tex`\\Theta_Q(L) = 1 + \\Theta_1 L^s \\cdots + \\Theta_Q L^{sQ}`} is the seasonal MA \nlag polynomial,
\n${tex`\\{ \\varepsilon_t \\}_{t=1}`} is an iid sequence from ${tex`N(0,\\sigma_\\varepsilon^2)`}.
\n
\n\nThe spectral density of a seasonal ARMA(P,Q) process is
\n\n${tex`\nf(\\omega)=\\frac{\\sigma^{2}}{2\\pi}\\left|\\frac{\\theta_q(e^{-i\\omega})}{\\phi_{p}(e^{-i\\omega})}\\right|^{2}\\left|\\frac{\\Theta_Q(e^{-i\\omega})}{\\Phi_{P}(e^{-i\\omega})}\\right|^{2}\\quad \\text{ for }\\omega\\in [-\\pi,\\pi)\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":463,"value":"#### Model settings","pinned":false,"hidden":false,"mode":"md"},{"id":162,"value":"viewof input = Inputs.form({\n s: Inputs.range([1, 365], {value: 12, step: 1, label: \"seasonal period\"}),\n p: Inputs.range([0, 4], {value: 1, step: 1, label: \"no. of AR\"}),\n q: Inputs.range([0, 4], {value: 0, step: 1, label: \"no. of MA\"}),\n P: Inputs.range([0, 4], {value: 1, step: 1, label: \"no. of seasonal AR\"}),\n Q: Inputs.range([0, 4], {value: 0, step: 1, label: \"no. of seasonal MA\"})\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":981,"value":"viewof n = Inputs.range([0, 1000], {value: 100, label: \"time series length\", step: 1})","pinned":false,"hidden":false,"mode":"ojs"},{"id":465,"value":"#### Parameter settings","pinned":false,"hidden":false,"mode":"md"},{"id":467,"value":"viewof σ = Inputs.range([0.001, 5], {value: 1, label: tex`\\text{error stdev, }\\sigma`, step: 0.01})","pinned":false,"hidden":false,"mode":"ojs"},{"id":282,"value":"$2a","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"$2b","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"$2c","pinned":false,"hidden":false,"mode":"ojs"},{"id":121,"value":"$2d","pinned":false,"hidden":false,"mode":"ojs"},{"id":781,"value":"viewof simulatebutton = Inputs.button(\"Simulate new data\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":935,"value":"viewof plotlog = Inputs.toggle({label: \"plot log f(ω)\", value: false})","pinned":false,"hidden":false,"mode":"ojs"},{"id":608,"value":"\n \n \n \n

\n ${specplot}\n

\n \n

\n ${pltacf}\n

\n \n","pinned":false,"hidden":false,"mode":"html"},{"id":881,"value":"

","pinned":false,"hidden":false,"mode":"html"},{"id":877,"value":"\n \n \n \n

\n \n ${timeplot}\n

\n \n","pinned":false,"hidden":false,"mode":"html"},{"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":617,"value":"#### Plots","pinned":false,"hidden":false,"mode":"md"},{"id":152,"value":"specplot = Plot.plot({\n y: {\n label: yaxislabel,\n //label: yaxislabel,\n transform: (x) => {if (plotlog){return Math.log(x)}else{return x}}\n //domain: [0,1.05*jstat.normal.pdf(mu, mu, sigma)]\n },\n x: {label: \"Angular frequency ω\"},\n title: title,\n marks: [\n Plot.ruleY([0]),\n Plot.ruleX([0]),\n Plot.lineY(specdensdata, {x: \"ω\", y: \"specdens\", stroke: \"orange\", strokeWidth: 2})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":493,"value":"timeplot = Plot.plot({\n title: md`Realized time series from ARMA(${input.p},${input.q})×(${input.P},${input.Q}) process`,\n x: {\n domain: [0,n]\n },\n marks: [\n Plot.ruleX([0]),\n Plot.ruleY([d3.min(simulatedData.map(d => d.x))]),\n Plot.lineY(simulatedData, {filter: d => d.time < n, x: \"time\", y: \"x\", stroke: \"steelblue\", tip: true})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":806,"value":"pltacf = Plot.plot({\n title: md`Estimated ACF of ARMA(${input.p},${input.q})×(${input.P},${input.Q}) process with season *s=*${s}`,\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\"}),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":968,"value":"import {computeACF} from \"@mattiasvillani/statistics-tools\"","pinned":true,"hidden":false,"mode":"ojs"},{"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":96,"value":"#### Data","pinned":false,"hidden":false,"mode":"md"},{"id":256,"value":"s = input.s","pinned":false,"hidden":false,"mode":"ojs"},{"id":648,"value":"p = input.p","pinned":false,"hidden":false,"mode":"ojs"},{"id":649,"value":"q = input.q","pinned":false,"hidden":false,"mode":"ojs"},{"id":644,"value":"P = input.P","pinned":false,"hidden":false,"mode":"ojs"},{"id":646,"value":"Q = input.Q","pinned":false,"hidden":false,"mode":"ojs"},{"id":307,"value":"ϕvect_value = {\n if (ϕvect[0] == \"No regular AR\"){\n return [0]\n }else{\n return ϕvect\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":391,"value":"θvect_value = {\n if (θvect[0] == \"No regular MA\"){\n return [0]\n }else{\n return θvect\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":392,"value":"Φvect_value = {\n if (Φvect[0] == \"No seasonal AR\"){\n return [0]\n }else{\n return Φvect\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":394,"value":"Θvect_value = {\n if (Θvect[0] == \"No seasonal MA\"){\n return [0]\n }else{\n return Θvect\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":99,"value":"specdensdata = {\n const ωgrid = d3.range(0,Math.PI, 0.001);\n return ωgrid.map(d => ({ω: d, specdens: SpecDens(d, ϕvect_value, θvect_value, Φvect_value, Θvect_value, σ, s)}))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":939,"value":"yaxislabel = {if (plotlog){return \"log Spectral density log f(ω)\"}else{return \"Spectral density f(ω)\"}}","pinned":false,"hidden":false,"mode":"ojs"},{"id":951,"value":"title = {if (plotlog){\n return md`Log spectral density of an ARMA(${input.p},${input.q})×(${input.P},${input.Q}) process with season *s=*${s}`\n }else{\n return md`Spectral density of an ARMA(${input.p},${input.q})×(${input.P},${input.Q}) process with season *s=*${s}`\n }\n }","pinned":false,"hidden":false,"mode":"ojs"},{"id":94,"value":"#### Functions","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"$2e","pinned":false,"hidden":false,"mode":"ojs"},{"id":514,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":685,"value":"function multiplyPolynomials(poly1, poly2) {\n const result = [];\n const degree1 = poly1.length - 1;\n const degree2 = poly2.length - 1;\n\n for (let i = 0; i <= degree1 + degree2; i++) {\n result[i] = 0;\n }\n\n for (let i = 0; i <= degree1; i++) {\n for (let j = 0; j <= degree2; j++) {\n result[i + j] += poly1[i] * poly2[j];\n }\n }\n\n return result;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":638,"value":"function combineARPoly(ϕvect, Φvect) {\n\n // No seasonal\n if (typeof Φvect[0] == \"string\"){\n if (ϕvect == \"No regular AR\"){\n return []\n }else{\n return ϕvect.map(d => -d)\n }\n }\n\n // No regular\n if (typeof ϕvect[0] == \"string\"){\n var combinedPoly = d3.range(P*s).map(d => 0)\n for (const [index, Φ] of Φvect.entries()){\n combinedPoly[s*(index + 1) - 1] = -Φ\n }\n return combinedPoly\n }\n\n // Both regular and seasonal\n var combinedPolyRegular = d3.range(1 + p).map(d => 0)\n var combinedPolySeasonal= d3.range(1 + P*s).map(d => 0)\n combinedPolyRegular[0] = 1\n combinedPolySeasonal[0] = 1\n for (const [index, ϕ] of ϕvect.entries()){\n combinedPolyRegular[1 + index] = -ϕ\n }\n for (const [index, Φ] of Φvect.entries()){\n combinedPolySeasonal[s*(index + 1)] = -Φ\n }\n \n return multiplyPolynomials(combinedPolyRegular, combinedPolySeasonal).slice(1)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":719,"value":"function combineMAPoly(θvect, Θvect) {\n\n // No seasonal\n if (typeof Θvect[0] == \"string\"){\n if (θvect == \"No regular MA\"){\n return []\n }else{\n return θvect\n }\n }\n\n // No regular\n if (typeof θvect[0] == \"string\"){\n var combinedPoly = d3.range(Q*s).map(d => 0)\n for (const [index, Θ] of Θvect.entries()){\n combinedPoly[s*(index + 1) - 1] = Θ\n }\n return combinedPoly\n }\n\n // Both regular and seasonal\n var combinedPolyRegular = d3.range(1 + q).map(d => 0)\n var combinedPolySeasonal= d3.range(1 + Q*s).map(d => 0)\n combinedPolyRegular[0] = 1\n combinedPolySeasonal[0] = 1\n for (const [index, θ] of θvect.entries()){\n combinedPolyRegular[1 + index] = θ\n }\n for (const [index, Θ] of Θvect.entries()){\n combinedPolySeasonal[s*(index + 1)] = Θ\n }\n \n return multiplyPolynomials(combinedPolyRegular, combinedPolySeasonal).slice(1)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":519,"value":"simulatedData = simARMA(1000, ϕvect, θvect, Φvect, Θvect, 0, σ)","pinned":false,"hidden":false,"mode":"ojs"},{"id":838,"value":"n_acf_lags = Math.max(3*(p + P*s + q + Q*s),10)","pinned":false,"hidden":false,"mode":"ojs"},{"id":804,"value":"acf = computeACF(simulatedData.map(d => d.x), n_acf_lags)","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W30","stdlib":"1","titling":"auto"}