Time Series:
Start = 1
End = 20
Frequency = 1
[1] -1.0615712 -3.3427201 -2.9095297 -1.5579658 -2.6942738 -3.6786334
[7] -3.3995961 -2.5009399 -3.6425151 -2.7769243 -2.9941547 -3.6164653
[13] -1.5792640 -1.6316549 0.2458082 -0.2222871 0.5333117 1.7562844
[19] 0.2204732 -0.3471451
The model
An autoregressive model of order \(p\), \(AR(p)\), is defined as:
\[{\displaystyle X_{t}=c+\sum _{i=1}^{p}\varphi _{i}X_{t-i}+\varepsilon _{t}\,}\]
Where \(\varepsilon_t\) is a White Noise, \(c\), \(\varphi_i\) parameters of the model. The structure of this model is easy of understand: The next value \(X_t\) is a linear combination of the past values plus a random noise -or innovation-.
In R we can simulate this model with stats::arima.sim
function for example, if we want simulate an \(AR(1)\) with \(\varphi _{1} = \varphi = 0.9\) with \(0\) mean:
\[{\displaystyle X_{t}= \varphi X_{t-i} + \varepsilon _{t} = 0.9 \times X_{t-i} + \varepsilon _{t}\,}\] We code:
Interactivity
There is a popular example in highcharts: Dynamic Update where each second a random value is added to the line series https://www.highcharts.com/demo/dynamic-update.
The key part of this example is the next function which every second generate a random value and then is added to the series.
Code
// This is Javascript code :) :/ :S XD
function () {
// set up the updating of the chart each second
var series = this.series[0];
setInterval(function () {
var x = (new Date()).getTime(), // current time
= Math.random();
y .addPoint([x, y], true, true);
series, 1000);
} }
We’ll take the example as a template and change to store the last values of the process and generate the random noise and the new value. To get the linear combination between the past values and the autoregressive parameter we used the dot product. For the generation of a normal random number we can use the the Box-Muller transform.
Let’s create some
Code
[1] -4.172 -6.993 -6.458 -6.460 -5.176 -5.944 -5.448 -5.063 -4.854 -5.728
The function would be something like this:
Code
load_fn <- "function () {{
dot = (a, b) => a.map((x, i) => a[i] * b[i]).reduce((m, n) => m + n);
function randn_bm() {{
var u = 0, v = 0;
while(u === 0) u = Math.random(); //Converting [0,1) to (0,1)
while(v === 0) v = Math.random();
return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v );
}}
var dat = [{ data }];
var ar = [{ ar }];
dat = dat.slice(-ar.length);
// set up the updating of the chart each second
var series = this.series[0];
setInterval(function () {{
console.log(dat);
var new_value = dot(dat, ar) + randn_bm();
new_value = Math.round(1000 * new_value)/1000
console.log(new_value)
dat.shift();
dat.push(new_value);
series.addPoint([new_value]);
//if (series.data.length < 500) {{
// series.addPoint([new_value], true, false);
//}} else {{
// series.addPoint([new_value], true, true);
//}}
}}, { time });
}}"
load_fn_glued <- str_glue(
load_fn,
data = str_c(ts, collapse = ","),
ar = str_c(ar, collapse = ","),
time = time * 1000
)
load_fn_glued
function () {
dot = (a, b) => a.map((x, i) => a[i] * b[i]).reduce((m, n) => m + n);
function randn_bm() {
var u = 0, v = 0;
while(u === 0) u = Math.random(); //Converting [0,1) to (0,1)
while(v === 0) v = Math.random();
return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v );
}
var dat = [-4.172,-6.993,-6.458,-6.46,-5.176,-5.944,-5.448,-5.063,-4.854,-5.728];
var ar = [0.85,-0.1,0.2];
dat = dat.slice(-ar.length);
// set up the updating of the chart each second
var series = this.series[0];
setInterval(function () {
console.log(dat);
var new_value = dot(dat, ar) + randn_bm();
new_value = Math.round(1000 * new_value)/1000
console.log(new_value)
dat.shift();
dat.push(new_value);
series.addPoint([new_value]);
//if (series.data.length < 500) {
// series.addPoint([new_value], true, false);
//} else {
// series.addPoint([new_value], true, true);
//}
}, 1000);
}
Code
And it works perfectly in the first try! Nah! This take me some time in the jsfiddle.
All according to keikaku.
Improving it
We can improve the chart doing for example:
Make a
plotLine
in y axis at the value 0 to show as reference.Make a second yAxis to the right showing the last generated value, this can be done using the next function in the
tickPositioner
argument.
Code
tick_post_fn <- "function(min,max){
var data = this.chart.yAxis[0].series[0].processedYData;
//last point
return [Math.round(1000 * data[data.length-1])/1000];
}"
Enable the navigator panel
hc_navigator(enabled = TRUE)
.We can put an informative title showing the specification of the simulated model.
Code
cat(formula)
\[ X_t = 0.85 \times X_{ t - 1 } -0.1 \times X_{ t - 2 } + 0.2 \times X_{ t - 3 } + \epsilon_t\]
Alternatively, formula using subscript and superscript tags.
Code
formula <- purrr::map2(ar, seq_along(ar), function(par, t){
if(par > 0 & t > 1) {
par <- str_c("+ ", par)
}
htmltools::tagList(par, "×", "X", tags$sub(stringr::str_glue("t - { i }", i = t)))
})
formula <- purrr::reduce(formula, htmltools::tagList)
formula <- str_c("X", tags$sub("t") %>% as.character(), " = ",formula %>% as.character(), " + ε", tags$sub("t") %>% as.character())
Code
cat(formula)
Xt = 0.85 × X t - 1 -0.1 × X t - 2 + 0.2 × X t - 3 + εt
- Helper button to remove last values.
Code
rm_poinst_fn <- "function () {
for (var i = 1; i <= 500; i++) {
if (this.series[0].data.length) {
this.series[0].data[0].remove();
}
}
}"
- And others tweaks.
Code
highchart() %>%
hc_add_series(data = ts, name = "Process") %>%
hc_chart(
events = list(load = JS(load_fn_glued)),
animation = list(duration = time*1000/2)
) %>%
hc_title(text = "Autoregressive process") %>%
hc_subtitle(text = as.character(formula), useHTML = TRUE) %>%
hc_plotOptions(series = list(marker = list(enabled = FALSE))) %>%
hc_tooltip(valueDecimals = 3) %>%
hc_xAxis(width = "95%") %>%
hc_exporting(
enabled = TRUE,
buttons = list(
list(
text = "Remove last 500 values",
onclick = JS(rm_poinst_fn),
theme = list(stroke = 'silver')
)
)
) %>%
hc_yAxis_multiples(
# default axis
list(
title = list(text = ""),
plotLines = list(
list(value = 0, width = 2, color = "#AAA", zIndex = 1)
)
),
# opposite axis
list(
title = list(text = ""),
linkedTo = 0,
opposite = TRUE,
gridLineWidth = 0,
tickPositioner = JS(tick_post_fn)
)
) %>%
hc_navigator(
enabled = TRUE,
series = list(type = "line"),
xAxis = list(labels = list(enabled = FALSE), width = "95%")
)
Extending it
We can encapsulate all this code into a function:
Code
library(highcharter)
library(stringr)
sim_ar_hc <- function(ar = c(0.3, 0.2), time = 1){
# the first values of the model
ts <- as.vector(arima.sim(model = list(ar = ar), n = 10))
ts <- round(ts, 3)
ts
load_fn <- "function () {{
dot = (a, b) => a.map((x, i) => a[i] * b[i]).reduce((m, n) => m + n);
function randn_bm() {{
var u = 0, v = 0;
while(u === 0) u = Math.random(); //Converting [0,1) to (0,1)
while(v === 0) v = Math.random();
return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v );
}}
var dat = [{ data }];
var ar = [{ ar }];
dat = dat.slice(-ar.length);
// set up the updating of the chart each second
var series = this.series[0];
setInterval(function () {{
console.log(dat);
var new_value = dot(dat, ar) + randn_bm();
new_value = Math.round(1000 * new_value)/1000
console.log(new_value)
dat.shift();
dat.push(new_value);
series.addPoint([new_value]);
//if (series.data.length < 500) {{
// series.addPoint([new_value], true, false);
//}} else {{
// series.addPoint([new_value], true, true);
//}}
}}, { time });
}}"
load_fn_glued <- str_glue(
load_fn,
data = str_c(ts, collapse = ","),
ar = str_c(ar, collapse = ","),
time = time * 1000
)
tick_post_fn <- "function(min,max){
var data = this.chart.yAxis[0].series[0].processedYData;
//last point
return [Math.round(1000 * data[data.length-1])/1000];
}"
rm_poinst_fn <- "function () {
for (var i = 1; i <= 500; i++) {
if (this.series[0].data.length) {
this.series[0].data[0].remove();
}
}
}"
formula <- purrr::map2(ar, seq_along(ar), function(par, t){
if(par > 0 & t > 1) {
par <- str_c("+ ", par)
}
htmltools::tagList(par, "×", "X", tags$sub(stringr::str_glue("t - { i }", i = t)))
})
formula <- purrr::reduce(formula, htmltools::tagList)
formula <- str_c("X", tags$sub("t") %>% as.character(), " = ",formula %>% as.character(), " + ε", tags$sub("t") %>% as.character())
hc <- highchart() %>%
hc_add_series(data = ts, name = "Process") %>%
hc_chart(
events = list(load = JS(load_fn_glued)),
animation = list(duration = time*1000/2)
) %>%
hc_title(text = "Autoregressive process") %>%
hc_subtitle(text = formula, useHTML = TRUE) %>%
hc_plotOptions(series = list(marker = list(enabled = FALSE))) %>%
hc_tooltip(valueDecimals = 3) %>%
hc_xAxis(width = "95%") %>%
hc_exporting(
enabled = TRUE,
buttons = list(
list(
text = "Remove last 500 values",
onclick = JS(rm_poinst_fn),
theme = list(stroke = 'silver')
)
)
) %>%
hc_yAxis_multiples(
# default axis
list(
title = list(text = ""),
plotLines = list(
list(value = 0, width = 2, color = "#AAA", zIndex = 1)
)
),
# opposite axis
list(
title = list(text = ""),
linkedTo = 0,
opposite = TRUE,
gridLineWidth = 0,
tickPositioner = JS(tick_post_fn)
)
) %>%
hc_navigator(
enabled = TRUE,
series = list(type = "line"),
xAxis = list(labels = list(enabled = FALSE), width = "95%")
)
hc
}
White Noise:
Code
sim_ar_hc(ar = c(0))
Some traditional AR model:
Code
sim_ar_hc(ar = c(.8, .1))
Let it shine (with {shiny} package)
In this case we can update the highcharts using the new set of proxy functions. You can check the code of the shiny app in https://github.com/jbkunst/shiny-apps-edu/ and the app live is in https://jbkunst.shinyapps.io/arma-process
References
- https://jsfiddle.net/gh/get/library/pure/highcharts/highcharts/tree/master/samples/highcharts/demo/dynamic-update
- https://jsfiddle.net/BlackLabel/evypfr1L/
- http://jsfiddle.net/upt4cbqj/
- http://www.java2s.com/Tutorials/highcharts/Example/Series_Data/Remove_first_data_item_from_series.htm
Reuse
Citation
@online{kunst_fuentes2021,
author = {Kunst Fuentes, Joshua},
title = {Simulating and Watching Autoregressive Processes},
date = {2021-06-18},
url = {https://jkunst.com/blog/posts/2021-07-17-simulating-and-watching-autoregressive-processes/},
langid = {en}
}