Data, Code and Visualization - Simulating and watching autoregressive processes

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

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();
        series.addPoint([x, y], true, true);
    }, 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] -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

highchart() %>% 
  hc_add_series(data = ts) %>% 
  hc_chart(
    events = list(load = JS(load_fn_glued)), 
    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)

X_t = 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(
      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

https://creativecommons.org/licenses/by/4.0/

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}
}

For attribution, please cite this work as:

Joshua Kunst Fuentes. 2021. “Simulating and Watching Autoregressive Processes.” June 18, 2021. https://jkunst.com/blog/posts/2021-07-17-simulating-and-watching-autoregressive-processes.