# Simulating and watching autoregressive processes

Always love time series when I studied statistics. This is a post to honor the time series (no {shiny} required).
minipost
stats
highcharts
highcharter
time-series
data-visualization
Author

Joshua Kunst Fuentes

Published

June 18, 2021

## 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:

Code
ar_model <-  list(ar = c(0.9))

arima.sim(model = ar_model, n = 20)
Time Series:
Start = 1
End = 20
Frequency = 1
[1] -1.77672932 -0.42137552 -1.14412458 -1.69266885 -0.89272259 -0.91866922
[7] -0.19412967  0.47567958 -0.26995749 -0.72483040  0.03042900 -0.31695851
[13]  1.04817286 -0.05238642 -0.49166004  0.01527420  0.52553654 -0.58810687
[19] -0.60426529 -0.75117670

## 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
y = Math.random();
}, 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
library(highcharter)
library(stringr)

# model
ar <-  c(0.85, -0.1, 0.2)

# time to update the chart in seconds
time <- 1

# the first values of the model
ts <- as.vector(arima.sim(model = list(ar = ar), n = 10))
ts <- round(ts, 3)
ts
 [1] 1.155 2.471 3.908 1.892 2.101 3.163 1.463 2.857 2.956 2.607

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);

//if (series.data.length < 500) {{
//}} else {{
//}}

}}, { time });

}}"

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 = [1.155,2.471,3.908,1.892,2.101,3.163,1.463,2.857,2.956,2.607];
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);

//if (series.data.length < 500) {
//} else {
//}

}, 1000);

}
Code
highchart() %>%
hc_chart(
animation = list(duration = time*1000/2)
)

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 tickPositionerargument.

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
formula <- purrr::map2(ar, seq_along(ar), function(par, t){

if(par > 0 & t > 1) {
par <- str_c("+ ", par)
}

stringr::str_glue("{ par } \\times X_{{ t - { i } }}", i = t)

})

formula <- purrr::reduce(formula, str_c, sep = " ")

formula <- str_c("$$X_t = ", formula, " + \\epsilon_t$$")
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(), " + &epsilon;",  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(), " + &epsilon;", tags$sub("t") %>% as.character())

hc <- highchart() %>%
hc_add_series(data = ts, name = "Process") %>%
hc_chart(
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 = ""),
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))

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

## Citation

BibTeX citation:
@online{kunstfuentes2021,
author = {Joshua Kunst Fuentes},
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}
}