Loading packages, defining colors and using data

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. If the variable export_graphs is set to TRUE, then the graphs will be exported as pdf-files. In addition, we define a set of colors here to make graphs look more beautiful. Furthermore, we define set the standard deviation of innovation shocks to an AR(1) process used in this section to a value of \(0.1\).

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

# should graphs be exported to pdf
export_pdf <- FALSE

# define some colors
mygreen <- "#00BA38"
myblue  <- "#619CFF"
myred   <- "#F8766D"
  
# standard deviation of epsilon
sd_eps = 0.1

Impulse response functions for different AR(1) coefficients

We now want to plot different impulse response functions of an AR(1) process with different degrees of persistence \(\rho\). To this end, we simulate the behavior of the AR(1) process for 20 periods. We therefore create an array x that contains each of these periods (quarters) and start a data frame dat that contains the periods in the first column. In a second step, we calculate the behavior of the AR(1) for persistence values in between \(-0.9\) and \(0.9\). We assume that the process receives a one standard deviation shocks (the impulse) at period \(t = 1\). In any successive period, there is no additional shock, but the original shock just fades out at the rate \(\rho\) (response).

# generate sequence of 20 quarters
x <- seq(1, 20, 1)
dat <- data.frame(quarter=x)

# iterate over different values
i <- 1
for(rho in c(-0.9, -0.5, 0, 0.5, 0.9)) {
  dat[paste("rho", i, sep="")] <- sd_eps*rho^c(0:(length(x)-1))
  i <- i+1
}


We now plot the impulse response functions for all non-negative values of \(\rho\) simulated above. It is easy to see that the persistence parameter \(\rho\) is the central variable that governs the memory of the process. High levels of persistence lead the initial shock to persist (to be remembered) for a high number of periods. Under a persistence of \(\rho = 0\), the shock only lasts for one period and then is immediately gone. Note that in the absence of shocks, the AR(1) process (with mean \(0\)) is just a first-order linear difference equation with \(a = \rho\) and \(b = 0\).

# create plot for positive rhos
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=quarter, y=rho3, color="l1"), linewidth=1) +
  geom_line(aes(x=quarter, y=rho4, color="l2"), linewidth=1) +
  geom_line(aes(x=quarter, y=rho5, color="l3"), linewidth=1) +
  scale_x_continuous(breaks=seq(0, 20, 4), expand=c(0, 0)) +  
  scale_y_continuous(breaks=seq(0, 0.1, 0.025)) +  
  labs(x = "Quarter t",
       y = expression(tilde(x)[t])) +
  scale_color_manual(breaks = c("l1", "l2", "l3"),
                     labels = c(expression(paste(rho," = 0.0")), expression(paste(rho," = 0.5")), expression(paste(rho," = 0.9"))),
                     values = c(mygreen, myblue, myred), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)


Remembering the convergence properties of the linear difference equation, it is not surprising that the impulse response function is alternating around zero when ever \(\rho < 0\). This can be seen from the following graph.

# create plot for negative rhos
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=quarter, y=rho2, color="l1"), linewidth=1) +
  geom_line(aes(x=quarter, y=rho1, color="l2"), linewidth=1) +
  scale_x_continuous(breaks=seq(0, 20, 4), expand=c(0, 0)) +  
  scale_y_continuous(breaks=seq(-0.1, 0.1, 0.025)) +  
  labs(x = "Quarter t",
       y = expression(tilde(x)[t])) +
  scale_color_manual(breaks = c("l1", "l2"),
                     labels = c(expression(paste(rho," = -0.5")), expression(paste(rho," = -0.9"))),
                     values = c(mygreen, myblue), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)

Distribution for \(\log(X_t)\) for different values of \(\rho\)

In this graph, we calculate the unconditional distribution (density function) for \(\log(X_t)\), which is just a normal distribution with mean \(\mu = 0\) and an unconditional variance of \(\sigma_2 = \frac{\sigma_\epsilon^2}{1-\rho^2}\), meaning a standard deviation of \(\sigma = \frac{\sigma_\epsilon}{\sqrt{1-\rho^2}}\). We can use the R function dnorm to calculate the density function of the normal distribution at different values x. We choose values of x ranging in between \(-1\) and \(1\) with a distance of \(0.01\). Finally, we plot the density functions for different values of \(\rho\). The larger is \(\rho\), the wider is the unconditional distribution of the process.

x = seq(-1, 1, 0.01)
dat <- data.frame(x)

# iterate over different values
i <- 1
for(rho in c(0, 0.5, 0.9)) {
  dat[paste("rho", i, sep="")] <- dnorm(x, mean=0, sd=sd_eps/sqrt(1-rho^2))
  i <- i+1
}

# create plot for different values of rho
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=x, y=rho1, color="l1"), linewidth=1) +
  geom_line(aes(x=x, y=rho2, color="l2"), linewidth=1) +
  geom_line(aes(x=x, y=rho3, color="l3"), linewidth=1) +
  scale_x_continuous(breaks=seq(-1, 1, 0.25), expand=c(0, 0)) +  
  labs(x = expression(log(x[t])),
       y = "Density") +
  scale_color_manual(breaks = c("l1", "l2", "l3"),
                     labels = c(expression(paste(rho," = 0.0")), expression(paste(rho," = 0.5")), expression(paste(rho," = 0.9"))),
                     values = c(mygreen, myblue, myred), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)

Distribution for \(X_t\) for different values of \(\rho\)

Now, we are interested in the distribution of \(X_t\) (not \(\log(X_t)\)). While \(\log(X_t)\) follows a simple normal distribution, \(X_t = \exp(\log(X_t))\) has to follow a log-normal distribution. The distribution can just be calculated from taking the exponential of x. The log-normal distribution is bound from below by the value \(0\). In addition, it is not symmetric anymore.

# create plot for different values of rho
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=exp(x), y=rho1, color="l1"), linewidth=1) +
  geom_line(aes(x=exp(x), y=rho2, color="l2"), linewidth=1) +
  geom_line(aes(x=exp(x), y=rho3, color="l3"), linewidth=1) +
  scale_x_continuous(breaks=seq(0, 3, 0.5), expand=c(0, 0)) +  
  labs(x = expression(x[t]),
       y = "Density") +
  scale_color_manual(breaks = c("l1", "l2", "l3"),
                     labels = c(expression(paste(rho," = 0.0")), expression(paste(rho," = 0.5")), expression(paste(rho," = 0.9"))),
                     values = c(mygreen, myblue, myred), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)

Distribution for \(\log(X_t)\) for different values of \(\bar X\)

Next, we keep the persistence parameter \(\rho = 0.9\) constant and change the long-run mean \(\bar X\) of the process. Increasing the mean shifts the distribution of \(\log(X_t)\) to the right.

rho = 0.9
x = seq(-1, 1.5, 0.01)
dat <- data.frame(x)

# iterate over different values
i <- 1
for(X_bar in c(0, 0.25, 0.5)) {
  y <- c(0:(length(x)-1))
  dat[paste("rho", i, sep="")] <- dnorm(x, mean=X_bar, sd=sd_eps/sqrt(1-rho^2))
  i <- i+1
}

# create plot for different values of rho
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=x, y=rho1, color="l1"), linewidth=1) +
  geom_line(aes(x=x, y=rho2, color="l2"), linewidth=1) +
  geom_line(aes(x=x, y=rho3, color="l3"), linewidth=1) +
  scale_x_continuous(breaks=seq(-1, 1.5, 0.25), expand=c(0, 0)) +  
  labs(x = expression(log(x[t])),
       y = "Density") +
  scale_color_manual(breaks = c("l1", "l2", "l3"),
                     labels = c(expression(paste(bar(X)," = 0.00")), expression(paste(bar(X)," = 0.25")), expression(paste(bar(X)," = 0.50"))),
                     values = c(mygreen, myblue, myred), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)

Distribution for \(X_t\) for different values of \(\bar X\)

The same is true for the distribution of \(X_t\).

# create plot for different values of rho
myplot <- ggplot(data = dat) + 
  geom_line(aes(x=exp(x), y=rho1, color="l1"), linewidth=1) +
  geom_line(aes(x=exp(x), y=rho2, color="l2"), linewidth=1) +
  geom_line(aes(x=exp(x), y=rho3, color="l3"), linewidth=1) +
  scale_x_continuous(breaks=seq(0, 4.5, 0.5), expand=c(0, 0)) +  
  
  labs(x = expression(x[t]),
       y = "Density") +
  scale_color_manual(breaks = c("l1", "l2", "l3"),
                     labels = c(expression(paste(bar(X)," = 0.00")), expression(paste(bar(X)," = 0.25")), expression(paste(bar(X)," = 0.50"))),
                     values = c(mygreen, myblue, myred), name='') +
  theme_bw() + 
  theme(legend.position="bottom")

# print the plot
print(myplot)