Loading packages etc.

We first clear the workspace using rm(list = ls())and then include all packages we need. If a package is missing in your R distribution (which is quite likely initially), just use install.packages("package_name") with the respective package name to install it on your system. If you execute the code in the file install_packages.R, then all necessary packages will be installed into your R distribution.

rm(list = ls())
library(reshape2)
library(base)
library(ggplot2)
library(grid)
library(stringr)
library(tidyverse)

A function for simulating a linear DE forward

Now we generate a function that allows us to simulate the DE forward. Input arguments to this function are a starting value y_0 as well as the parameters of the difference equation a and b. Finally, the function should simulate T periods of the difference equation. The function works as follows:

  1. We set the first value of an array y to y_0.

  2. We then iterate over all dates in time t and calculate the next value y[t] from the previous value using the specific form of the linear DE. Note that the standard index of a sequence in R starts at 1 and not at 0. Therefore we have to calculate T+1 values if we want to have the sequences running from 0 to T.

Finally, we return the sequence y as a result of the function.

diff_eq <- function(y_0, a, b, T) {
  y <- y_0
  for (t in 2:(T+1)) {
      y[t] = a*y[t-1] + b
  }
  return(y)
}

Plotting the solution to the DE

Now we can generate a function that allows us to plot the resulting forward simulated DE. The function receives the same arguments as the function diff_eq shown above. It directly makes use of this function and simulates the solution to the difference equation. We then store the resulting data as well as an array that contains the relevant time dimension it a data frame called dat. We then use ggplot2 to generate a nice-looking plot. This requires a bit of skill. But there are ample ggplot tutorials out there that can help, like this one.

plot_diff <- function(y_0, a, b, T) {
  
  # simulate equation forward
  y <- diff_eq(y_0, a, b, T)
  
  # generate plotting data
  dat <- data.frame(
    time = seq(0, T, 1),
    DE = y
  )
  
  # now generate the plot
  myplot <- ggplot(data=dat) + 
    geom_line(aes(x=time, y=DE), color="darkblue", linewidth=1) +
    geom_point(aes(x=time, y=DE), color="darkblue", size=3, show.legend = FALSE) +
    coord_cartesian(xlim=c(0, T)) + 
    scale_x_continuous(breaks=seq(0, T, 2)) +
    labs(x = "Time t",
         y = expression(paste("Solution of linear DE ", y[t])),
         title= bquote("Linear Difference Equation: " ~ a == .(a) ~ " , " ~ b == .(b) ~ " , " ~ y[0] == .(y_0))) +
    theme_bw()
  
  # print the plot
  print(myplot)
}

Some Examples

Now we can start generating plots. We want to simulate 15 periods of the solution to the linear difference equation. We therefore define a variable T = 14. We then define different parameter values for a, b and y_0 and want to see how the solution to the DE behaves. Setting a = 2 b = 10 and y_0 = 0, we obviously obtain a solution that diverges.

T = 14;

a = 2 ; b = 10 ; y_0 = 0
plot_diff(y_0, a, b, T)


Changing a to a value of -2 will lead the solution to diverge again, but in an alternating way.

a = -2 ; b = 10 ; y_0 = 0
plot_diff(y_0, a, b, T)


With a value of a = 1 the solution of the difference equation evolves linearly over time with a slope of b, i.e. we have \(y_t = y_0 + b\cdot t\).

a = 1; b = 10 ; y_0 = 0
plot_diff(y_0, a, b, T)


With a = -1, the sequence alternates between the two values \(0\) and \(10\).

a = -1; b = 10 ; y_0 = 0
plot_diff(y_0, a, b, T)

Adding the steady state

In a next step, we want to study convergence properties relative to the steady state. We therefore want to augment our plot to function by adding the steady state value \(\bar y = \frac{b}{1-a}\) to the plot. The rest of the function is as above.

plot_diff_steady <- function(y_0, a, b, T) {

  # simulate equation forward
  y <- diff_eq(y_0, a, b, T)
  
  # generate plotting data
  dat <- data.frame(
    time = seq(0, T, 1),
    DE = y,
    stst= array(b/(1-a), c(T+1))
  )
  
  # now generate the plot
  myplot <- ggplot(data=dat) + 
    geom_line(aes(x=time, y=stst, color="Steady State"), linewidth=1) +
    geom_line(aes(x=time, y=DE, color="Solution DE"), linewidth=1) +
    geom_point(aes(x=time, y=DE), color="darkblue", size=3, show.legend = FALSE) +
    scale_color_manual(name = "", values = c("Steady State"="red", "Solution DE"="darkblue")) +
    coord_cartesian(xlim=c(0, T)) + 
    scale_x_continuous(breaks=seq(0, T, 2)) +
    labs(x = "Time t",
         y = expression(paste("Solution of linear DE ", y[t])),
         title= bquote("Linear Difference Equation: " ~ a == .(a) ~ " , " ~ b == .(b) ~ " , " ~ y[0] == .(y_0))) +
    theme_bw()
  
  # print the plot
  print(myplot)
}


Now let us see again what happens for a = 2. The sequence again diverges.

a = 2 ; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)


But if we start directly in the steady state y_0 = b/(1-a) = -10, the system is in rest.

a = 2 ; b = 10 ; y_0 = b/(1-a)
plot_diff_steady(y_0, a, b, T)


Yet, if we depart from the steady state only a tiny little bit, we again get a diverging sequence.

a = 2 ; b = 10 ; y_0 = b/(1-a)+0.01
plot_diff_steady(y_0, a, b, T)


The same is true for a = -1’. Here the system rests if we start in the steady state value of 5.

a = -1; b = 10 ; y_0 = b/(1-a)
plot_diff_steady(y_0, a, b, T)


But immediately alternates around the steady state once we depart a little bit from the stady state value.

a = -1; b = 10 ; y_0 = b/(1-a)-0.01
plot_diff_steady(y_0, a, b, T)


Now if we choose an a with \(|a| < 1\), we get a converging solution.

a = 0.75; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)


Convergence is faster, the smaller is a.

a = 0.50; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)

a = 0.25; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)


The same is true for values of \(a \in (-1, 0)\). But convergence to the steady state is alternating.

a = -0.75; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)

a = -0.50; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)

a = -0.25; b = 10 ; y_0 = 0
plot_diff_steady(y_0, a, b, T)