library(dplyr)
library(ggplot2)

Inspired by Tavish Srivastava’s step by step guide to learn time series modelling. \(x_t\) follows a data generating process based on the normal distribution. The data generating process depends on the previous observation plus an error term: \[x_t = \rho x_{t-1} + e_t\] When \(\rho=0\), \(x_t\) is a stationary series, drawn from a normal distribution. When \(\rho=1\), \(x_t\) is a random walk. We are interested in what is happening when \(\rho\) varies.

By adding a constant parameter, time series can also be modeled as such: \[x_t = \beta_0 + \beta_1 x_{t-1} + e_t\]

Processes where \(\beta_1=1\) are said to have a unit root. Below we kept \(\beta_0=0\) and used the \(\rho\) notation instead of \(\beta_1\).

Stationary series and random walk

rho=0 use rnorm() to generate a stationary series

er <- rnorm(1:300)
x <- er
plot(er, type = "l")

rho=1 use rnorm() and cumsum() to generate a random walk

Repeat the random number generation 4 times to see different version of the plot.

plot(cumsum(rnorm(1:300)), type="l")

plot(cumsum(rnorm(1:300)), type="l")

plot(cumsum(rnorm(1:300)), type="l")

plot(cumsum(rnorm(1:300)), type="l")

Introduce the rho parameter

Compare 6 values of rho, use a loop to generate \(x_t\)

Use the same error series as above.

autoregressive <- function(rho, error){
  x <- 0
  for (i in c(1:300)){
    x <- c(x, rho * x[i] + error[i])
    }
  return(x)
}

plot(autoregressive(rho=0, error=er), type ="l")

plot(autoregressive(rho=0.5, error=er), type ="l")

plot(autoregressive(rho=0.8, error=er), type ="l")

plot(autoregressive(rho=0.9, error=er), type ="l")

plot(autoregressive(rho=0.99, error=er), type ="l")

plot(autoregressive(rho=1, error=er), type ="l")

rho>0 use Reduce to generate \(x_t\) on a vector

Could I use Reduce() ? Examples from the Reduce documentation

## A simple function to compute continued fractions:
cfrac <- function(x) Reduce(function(u, v) u + 1 / v, x, right = TRUE)
## Continued fraction approximation for pi:
cfrac(c(3, 7, 15, 1, 292))
## A general-purpose cumulative adder:
cadd <- function(x) Reduce("+", x, accumulate = TRUE)
cadd(seq_len(7))

Use Reduce(accumulate=TRUE) for the data generating process (DGP) on a vector.

dgp <- function(x, rho) Reduce(function(u, v) rho * u + v, x, accumulate = TRUE)
x1 <- rnorm(1000)
plot(dgp(rnorm(x1), rho = 1), type = "l")
lines(x1)

# Use DGP on a data frame
dtf <- data.frame(a = rnorm(1e3)) %>%
  mutate(b = dgp(a, rho = 0.9))
plot(dtf$b, type = "l")

plot(dtf$a, dtf$b, type = "l")

Compare 6 values of rho, use Reduce to generate \(x_t\) on a data frame

Show a visualisation of a random error and the data generating process for different values of rho.

### Visualise several versions of the rho paramater in one data frame
#' This function will be used as such: group_by(rho) %>% do(dgp2()).
#' Therefore dtf should contain only one rho parameter.
#' @param dtf a data frame containing one rho parameter, 
dgp2 <- function(dtf){ 
  rho <- unique(dtf$rho)
  stopifnot(length(rho)==1)
  dtf %>%
    mutate(xt = Reduce(function(u, v) rho * u + v, x, accumulate = TRUE))
}

randomdata <- function(n=1000){
  merge(data.frame(x = rnorm(n),
                   rown = c(1:1000)),
        data.frame(n = 100,
                   rho = c(0,0.5,0.9,0.99,1,1.01))) %>%
    group_by(rho) %>%
    do(dgp2(.)) %>%
    data.frame %>%
    mutate(rho = paste("rho =",rho))
  }

ggplot(randomdata(), aes(x = rown, y = xt)) + 
  geom_line() + facet_wrap(~rho, scales = "free_y")

ggplot(randomdata(), aes(x = rown, y = xt)) + 
  geom_line() + facet_wrap(~rho, scales = "free_y")

You could plot different simulation on one plot with different colors.

randomdata2 <- rbind(mutate(randomdata(), run = 1),
                     mutate(randomdata(), run = 2),
                     mutate(randomdata(), run = 3),
                     mutate(randomdata(), run = 4))
ggplot(randomdata2, aes(x = rown, y = xt, color=as.factor(run))) + 
  geom_line() + facet_wrap(~rho, scales = "free_y")

Autocorrelation plot vor various values of rho

These plot are based on the same randomdata2 dataset generated for the plot above. Use the acf function. Use of ggplot2::geom_segment to generate acf plots was inspired by this question on Stackoverflow.

#' Calculate the autocorrelation estimates for a vector in a data frame
#' To be used after group_by(something) %>%  do(acf2(.)) 
#' @param dtf a dataframe containing xt data 
acf2 <- function(dtf){
  autocorr <- acf(dtf$xt, plot = FALSE)
  data.frame(acf = autocorr$acf,
             lag = autocorr$lag)
}
randomdata2 %>% 
  group_by(run, rho) %>%
  do(acf2(.)) %>%
  ggplot(aes(x = lag, y = acf)) +
  geom_hline(aes(yintercept = 0)) +
  geom_segment(aes(xend = lag, yend = 0) ) + 
  facet_grid(run ~ rho)