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)
library(gridExtra)

A function for simulating a non-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 curvature parameter of the function that governs the non-linear difference equation \(f(y) = y^\alpha\). Finally, the function should simulate T periods of the difference equation.

diff_eq <- function(y_0, alpha, T) {
  y <- y_0
  if(T > 0) {
    for (t in 2:(T+1)) {
      y[t] = y[t-1]^alpha
    }
  }
  return(y)
}

Generating a phase diagram

We next provide a function plot_phase that allows to draw a phase diagram of the difference equation specified in the above function. It receives the same input arguments as the function above. This function proceeds in several steps to construct the phase diagram:

  1. In a first step, we simulate the difference equation forward until date T.

  2. Next, we want to plot the functional form of the difference equation as well as a \(45°\)-line. We therefore construct a data frame that contains an equally spaced set of points on the interval \([0, 2]\) as well as the corresponding function values of \(f(y) = y^\alpha\). We plot the function and the \(45°\)-line by means of ggplot.

  3. Next we successively add the arrows for each step of the phase diagram. We do this using a for-loop that allows us to repeat the same steps over and over again. For each step of the phase diagram, we have to draw a vertical line (the arrow from the \(45°\)-line to the function \(f(y)\)) as well as a horizontal line (the arrow from the the function \(f(y)\) back to the \(45°\)-line). For each of the arrows, we construct a separate data frame dat and use the geom_path plot statement to add the arrow to the plot.

Note that, when all T arrow sets are added to the plot, we do not print the plot, but return the plot to the function that will call plot_phase, as we want to combine a phase diagram as well as the time diagram of the difference equation into one plot.

plot_phase <- function(y_0, alpha, T) {
  
  # simulate equation forward
  y <- diff_eq(y_0, alpha, T)
  
  # First plot the function itself
  dat <- data.frame(
    x = seq(0.0, 2, by=0.01),
    f = seq(0.0, 2, by=0.01)^alpha
  )
  
  myplot <- ggplot(data=dat) + 
    geom_line(aes(x=x, y=f), color="darkblue", linewidth=1) +
    geom_line(aes(x=x, y=x), color="darkblue", linewidth=0.5) +
    coord_cartesian(xlim=c(0, 2), ylim=c(0, 2)) + 
    scale_x_continuous(breaks=seq(0, 2, 0.2)) +
    labs(x = expression('y'['t']),
         y = expression('y'['t+1']),
         title= bquote("Phase Diagram: " ~ alpha == .(alpha) ~ " , " ~ y[0] == .(y_0))) +
    theme_bw()
  
  # add single steps to the phase diagram
  if(T >= 1) {
    for (t in 1:T) {
      
      # add the vertical line
      dat <- data.frame(d1 = c(y[t], y[t]),
                        d2 = c(y[t], y[t+1]))
      myplot <- myplot + geom_path(data=dat, aes(x=d1, y=d2), color = "red", linewidth=0.5, arrow = arrow(length = unit(0.02, "npc")))
      
      # add the horizontal line
      dat <- data.frame(d1 = c(y[t], y[t+1]),
                        d2 = c(y[t+1], y[t+1]))
      myplot <- myplot + geom_path(data=dat, aes(x=d1, y=d2), color = "red", linewidth=0.5, arrow = arrow(length = unit(0.02, "npc")))
    }
  }
  
  # print the plot
  return(myplot)
}

Plotting the time diagram of the solution to the DE

In addition to the phase diagram, we again want to plot a time diagram, as we did for linear difference equations. The function plot_diff is almost identical to the one we used for plotting linear DEs. The only differences are in the input values as well as the fact that we do not print the plot, but return the created plot through the return statement.

plot_diff <- function(y_0, alpha, T) {
  
  # simulate equation forward
  y <- diff_eq(y_0, alpha, T)
  
  # generate plotting data
  dat <- data.frame(
    time = seq(0, T, 1),
    DE = y,
    stst = array(1, c(T+1))
  )
  
  aval <- sprintf("%5.2f", alpha)
  yval <- sprintf("%5.2f", y_0)
  
  # now generate the plot
  myplot <- ggplot(data=dat) + 
    geom_line(aes(x=time, y=stst), color="red", linewidth=1) +
    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), ylim=c(0,2)) + 
    scale_x_continuous(breaks=seq(0, T, 2)) +
    labs(x = "Time t",
         y = expression(paste("Solution of non-linear DE ", y[t])),
         title= bquote("Non-Linear Difference Equation: " ~ alpha == .(alpha) ~ " , " ~ y[0] == .(y_0))) +
    theme_bw()
  
  # print the plot
  return(myplot)
}

Printing two plots next to each other

Now we want to print two plots next to each other. To do this, we first have to generate the phase diagram as well as the time diagram using the functions plot_phase and time_diagram we constructed above. To show two graphs next to each other, we can use the function grid.arrange out of the gridExtra package. This function takes the plots we want to show. The statement ncol=2 ensures that the two plots are arranged in two different columns, meaning next to each other.

plot_both <- function(y_0, alpha, T) {
  
  #  phase diagram
  phase_diagram <- plot_phase(y_0, alpha, T)
  
  # time diagram
  time_diagram <- plot_diff(y_0, alpha, T)
  
  # plot one next to the other
  grid.arrange(phase_diagram, time_diagram, ncol=2)
}

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 = 15. We then define different parameter values for alpha and y_0 and want to see how the solution to the DE behaves. We start out with \(\alpha = 0.5\) and a starting value that ensures monotone convergence to the steady state from the left.

T = 14;

plot_both(0.25, 0.5, T)


When we chose a starting point on the right of the steady state of \(\bar y = 1\), we get monotone convergence from the right.

plot_both(1.75, 0.5, T)


With a value of \(\alpha \in (-1,0)\), we get alternating convergence to the steady state.

plot_both(0.25, -0.5, T)


For \(\alpha > 1\), the steady state at \(\bar y = 1\) is instable. But the steady state at \(\bar y = 0\) is locally stable for starting value \(y_0 \in (0, 1)\).

plot_both(0.9, 1.5, T)


If we start above the steady state of \(\bar y = 1\), the sequence diverges monotonically.

plot_both(1.01, 1.5, T)


Finally, the sequence diverges in an alternating way for \(\alpha < -1\).

plot_both(1.01, -1.5, T)