Goals

Stationarity makes progress

For some reason it occurred to me that a lot of the analyses I was conducting were spurious, not because of the data, but because of my ineptitude at recognising a glaring issue with the data wrangling process… It wasn’t stationary lol. Hopefully we can figure that out in this attempt. The goal is to use data on the South African economy via the samadb package, to paint a picture of the relationships between different macroeconomic variables. We’ll be constructing our own monetary index, of the törnqvist variety, and looking into the relationship between that, the rate of unemployment, the rate of inflation and economic growth. We hope to be more purposeful in our approach to VAR models, investigating methods of de-trending data to gain more significant results before extending the analysis to a structural vector auto-regression. We start by clearing our global environment and loading our suite of packages (we probably won’t use all of them).

#clear our environment and any plots
rm(list = ls())
graphics.off()


# install.packages("xts")
library(xts)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# install.packages("PerformanceAnalytics")
library(PerformanceAnalytics)

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

# install.packages("vars")
library(vars)

## Loading required package: MASS

## Loading required package: strucchange

## Loading required package: sandwich

## Loading required package: urca

## Loading required package: lmtest

# install.packages("dplyr")
library(dplyr)

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:xts':
## 
##     first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#install.packages("samadb")
library(samadb)

## samadb 0.3.1, see help(samadb). Commercial users require a license at https://econdata.co.za/support/

# install.packages("TSA")
library(TSA)

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

# install.packages("stats")
library(stats)
# install.packages("ggplot2")
library(ggplot2)
# install.packages("imputeTS")
library(imputeTS)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA

## 
## Attaching package: 'imputeTS'

## The following object is masked from 'package:zoo':
## 
##     na.locf

# install.packages("lmom")
library(lmom)
# install.packages("zoo")
library(zoo)
# install.packages("tidyverse")
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ stringr::boundary() masks strucchange::boundary()
## ✖ dplyr::filter()     masks stats::filter()
## ✖ dplyr::first()      masks xts::first()
## ✖ dplyr::lag()        masks stats::lag()
## ✖ dplyr::last()       masks xts::last()
## ✖ dplyr::select()     masks MASS::select()
## ✖ readr::spec()       masks TSA::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# install.packages("matrixStats")
library(matrixStats)

## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:dplyr':
## 
##     count

# install.packages("tidyquant")
library(tidyquant)

## ── Attaching core tidyquant packages ─────────────────────── tidyquant 1.0.11 ──
## ✔ quantmod 0.4.26     ✔ TTR      0.24.4
## ── Conflicts ────────────────────────────────────────── tidyquant_conflicts() ──
## ✖ zoo::as.Date()                 masks base::as.Date()
## ✖ zoo::as.Date.numeric()         masks base::as.Date.numeric()
## ✖ stringr::boundary()            masks strucchange::boundary()
## ✖ matrixStats::count()           masks dplyr::count()
## ✖ dplyr::filter()                masks stats::filter()
## ✖ dplyr::first()                 masks xts::first()
## ✖ TSA::kurtosis()                masks PerformanceAnalytics::kurtosis()
## ✖ dplyr::lag()                   masks stats::lag()
## ✖ dplyr::last()                  masks xts::last()
## ✖ PerformanceAnalytics::legend() masks graphics::legend()
## ✖ imputeTS::na.locf()            masks zoo::na.locf()
## ✖ dplyr::select()                masks MASS::select()
## ✖ TSA::skewness()                masks PerformanceAnalytics::skewness()
## ✖ readr::spec()                  masks TSA::spec()
## ✖ quantmod::summary()            masks urca::summary(), base::summary()
## ✖ tidyquant::VAR()               masks vars::VAR()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## 
## Attaching package: 'tidyquant'
## 
## 
## The following object is masked from 'package:vars':
## 
##     VAR

# install.packages("IndexNumR")
library(IndexNumR)
#install.packages("tempdisagg")
library(tempdisagg)
#install.packages("urca")
library(urca)
#install.packages("lubridate")
library(lubridate)
#install.packages("seasonal)
library(seasonal)

## Warning in system(paste(x13.bin, file.path(tdir, "Testairline")), intern =
## TRUE, : running command
## '/Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/x13binary/bin/x13ashtml
## /var/folders/pn/jl67ndws11b3dw5mstg6xx6r0000gn/T//Rtmp2iOrrD/x13binary__d0fd7a1ba102__dir/Testairline
## 2>/dev/null' had status 134

## The binaries provided by 'x13binary' do not work on this
## machine. To get more information, run:
##   x13binary::checkX13binary()
## 
## You can set 'X13_PATH' manually if you intend to use your own
## binaries. See ?seasonal for details.
## 
## 
## Attaching package: 'seasonal'
## 
## The following object is masked from 'package:tibble':
## 
##     view

#install.packages("zoom")
library(zoom)

 library(tsbox)

Samadb and retrieving data

The specific functions that are important for our purposes are the sm_series() and sm_dataset() functions. First, we retrieve a dataframe of all the time series data, their number of observations and the specific “series” it’s retrieved from (“QB” for Quarterly Bulletin, provided by SARS, will be used primarily for the monetary aggregates).

df.allseries <- sm_series() ## there are 9526 data series in all. some of which have data going back into the 1960s.

Using this as our guide, we then go on to retrieve all the data-sets we need from the database. This part will look a lot like the previous iterations, with improvements to the naming. We’re using data from 2013 through all of 2019 in order to avoid the covid period, which will complicate our processes, and include data on own rates of return for many of our types of money. This data seems to have only begun being collected in 2013, or at least that is the earliest I could find meaningful observations (the rest of the observations were 1s which intuitively doesn’t seem like reliable estimates of the varying rates of interest.

##################################################################################
######Quarterly QB series

df.QuartQB <- samadb::sm_data(dsid = "QB",
                         series = "KBP1360M",
                         from = "2013-01-01",  # Dates must be strings
                         to = "2023-01-01",    # Dates must be strings
                         freq = "Q",
                         labels = TRUE,         # Use logicals
                         wide = TRUE,
                         ordered = TRUE)

##################################################################################
###############Quarterly National Accounts  dataseries
df.QBweekly <- samadb::sm_data(dsid = "QB",
                        series =  "KBP1405W" ,
                        from = "2013-01-01",  # Dates must be strings
                        to = "2023-01-01",    # Dates must be strings
                        freq = "W",
                        labels = TRUE,         # Use logicals
                        wide = TRUE,
                        ordered = TRUE)
#################################################################################
################ Monthly QB dataseries

df.QBMonbank <- samadb::sm_data(dsid = "QB_MONBANK",
                         series = "KBP1486M_M_N_N_LVL",
                         from = "2013-01-01",  # Dates must be strings
                         to = "2023-01-01",    # Dates must be strings
                         freq = "M",
                         labels = TRUE,         # Use logicals
                         wide = TRUE,
                         ordered = TRUE) ## Returns all monthly series in QB_MONBANK
################################################################################


df.MonthlyQB <- samadb::sm_data(dsid = "QB",
                         series = "KBP1360M",
                         from = "2013-01-01",  # Dates must be strings
                         to = "2023-01-01",    # Dates must be strings
                         freq = "M",
                         labels = TRUE,         # Use logicals
                         wide = TRUE,
                         ordered = TRUE) ## Returns all monthly series in QB
########################################################################################################

As before, we then select our inputs for the monetary agggregate:

df.MoneySupply<- df.MonthlyQB%>%
                  dplyr::select(date,
                                      KBP1312M, ## 1 CASH AND COIN CIRCULATING (own rate of return is negative CPI inflation/12)
                                      KBP1313M, ## 2 Cheque and Transmission Deposits
                                      KBP1314M, ## 3 DEMAND DEPOSITS (orr is KBP1601M)
                                      KBP1321M, ## 4 SAVINGS DEPOSITS (orr KBP1606M)
                                      KBP1316M, ## 5 SHORT-TERM DEPOSITS (orr is KBP1609M)
                                      KBP1322M, ## 6 MEDIUM-TERM DEPOSITS (orr is KBP2007M / KBP1611M)
                                      KBP1319M, ## 7 LONG-TERM DEPOSITS (orr is KBP1612M)
                                      
                                      KBP1070M,  ## 8 Cash managed, cheque and transmission deposits
                                      KBP1082M, ## 9 Debt securities
                                      
                                
                                      
                                      KBP1369M, ## 10 LOANS and ADVANCES
                                      
                                      KBP1334M, ## 11 Domestic capital and reserves 
                                      
                                      KBP1155M, ## - 12 DEPOSITS BY MONETARY INSTITUTIONS
                                      KBP1161M, ## 13 DISCOUNTED TREASURY BILLS
                                      
                                      KBP1200M, ## 14 MUTUAL BANK TRANSMISSION DEPOSITS
                                      KBP1201M, ## 15 M BANK SAVINGS DEPOSITS
                                      KBP1202M, ## 16 M BANK OTHER SHORT-MED TERM DEPOSITS
                                      KBP1203M, ## 17 M BANK LONG TERM DEPOSITS
                                      
                                      KBP1209M, ## 18 POSTBANK DEPOSITS
                                      )

For the sake of our internal rate of return, the first column has to be the negative of the rate of inflation. As such we begin by extracting a dataset for inflation seperately to the other internal rates of return so we can look at inflation on it’s own.

df.pi <- df.MonthlyQB %>%
  dplyr::select(
    date,
    KBP7146A # 1 CPI
  )

We can then go on to retrieve the rest of our rates of returns:

df.int <- df.QBMonbank %>%
  dplyr::select(
    
    KBP1604M_M_N_N_LVL, ## 2 BA930 Deposit Rates: Household Sector: Cheque Accounts
    KBP1605M_M_N_N_LVL, ## 3 BA930 Deposit Rates: Household Sector: Call Deposits 
    KBP1612M_M_N_N_LVL, ## 4 BA930 Deposit Rates: Household Sector: Fixed Deposits(5yrs+) Weighted Average (%). , Nominal
    KBP1609M_M_N_N_LVL, ## 5 BA930 Deposit Rates: Household Sector: Fixed Deposits(0-1yr)
    KBP1610M_M_N_N_LVL, ## 6 BA930 Deposit Rates: Household Sector: Fixed Deposits (1<x<3)
    KBP1611M_M_N_N_LVL, ## 7 BA930 Deposit Rates: Household Sector: Fixed Deposits(3yrs+>&<5yrs) Weighted Average (%). , Nominal
    KBP1616M_M_N_N_LVL, ## 8 BA930 Deposit Rates: Corporate Sector: Call Deposits Weighted Average (%). , Nominal
    KBP1603M_M_N_N_LVL, ## 9 Yield on Bonds Traded on the Stock Exchange: Government Bonds (20-30yrs)
    KBP1626M_M_N_N_LVL, ## 10 BA930 Lending Rates: All Domestic Private Sectors: Leasing Transactions: Flexible Rate Weighted Average. , Nominal
    KBP1606M_M_N_N_LVL, ## 11 BA930 Deposit Rates: Household Sector: Notice Deposit
    KBP1634M_M_N_N_LVL, ## 12 BA930 Lending Rates: Household Sector: Instalment Sale 

    KBP1618M_M_N_N_LVL, ## 13 Prime Overdraft Rate. End of Period, Nominal, 
    
    KBP1615M_M_N_N_LVL, ## 14 BA930 Deposit Rates: Corporate Sector: Cheque Accounts Weighted Average (%). , Nominal

    KBP1617M_M_N_N_LVL, ## 15 BA930 Deposit Rates: Corporate Sector: Notice Deposits (1-32 days)
    KBP1620M_M_N_N_LVL, ## 16 BA930 Deposit Rates: Corporate Sector: Fixed Deposits
    KBP1622M_M_N_N_LVL, ## 17 BA930 Deposit Rates: Corporate Sector: Fixed Deposits (3yrs & More, but less than 5yrs) Weighted Average (%). , Nominal
    KBP1602M_M_N_N_LVL ## 18 BA930 Deposit Rates: All Domestic Private Sectors: Other Weighted Average (%). , Nominal

  )

monthly_xts <- xts(, order.by = df.MoneySupply$date)## an empty xts obtect with the monthly index used in the money supply estimate. We'll use this to dissagregate quarterly data into monthly data as well as for our indexing purposes.

We can then combine our incomplete internal rates of return data with the inflation data:

df.irr <- cbind(df.pi[,-1],df.int)

We should probably collect the rest of our data at some point and TNTLTP so here we go.

Pause, real variables.

we want data on unemployment, nominal GDP, and real GDP.

df.unem <- df.QuartQB%>%
  dplyr::select(date,
                KBP7019K)
#////////////////////////////////////////////////////////////
df.NGDP <- df.QuartQB%>%
  dplyr::select(date,
                KBP6006L)
#////////////////////////////////////////////////////////////
df.RGDP <- df.QuartQB%>%
  dplyr::select(date,
                KBP6006D)

unemployment.xts <- xts(df.unem[,-1], order.by = df.unem$date)

nominalGDP.xts <- xts(df.NGDP[,-1], order.by = df.NGDP$date)

realGDP.xts <- xts(df.RGDP[,-1], order.by = df.RGDP$date)

This data is indexed quarterly but we want monthly observations. This means we’ll have to disaggregate our datatset.

First we combine our monthly data with our empty xts object, this will create a lot of blank entries (NAs) inbetween the quarterly observations.

While crude, we’ll use the na_interpolation function to approximate the intermediate values.

realGDP.xts <- cbind(monthly_xts, realGDP.xts)%>%
  na_kalman()

unemployment.xts <- cbind(monthly_xts, unemployment.xts)%>%
  na_kalman()

nominalGDP.xts <- cbind(monthly_xts, nominalGDP.xts)%>%
  na_kalman()

While we’re at it we might as well create time series objects of our monthly data.

MonQuant.xts <- xts(df.MoneySupply[,-1], order.by = df.MoneySupply$date) ## turn em to an xts object

MonQuant.xts[, 12] <- MonQuant.xts[, 12] * -1 ## Deposits by other monetary institutions is a negative to the total money supply

irr.xts <- xts(df.irr, order.by = df.pi$date)

irr.xts$KBP7146A <- irr.xts$KBP7146A *-(1/12)

GDPdeflator.xts <- (nominalGDP.xts/realGDP.xts)*100

Our log-transformed data looks like this:

unemployment.xts <- log(unemployment.xts)

realGDP.xts <- log(realGDP.xts)

nominalGDP.xts <- log(nominalGDP.xts)

GDPdeflator.xts <- log(GDPdeflator.xts)

pi.xts <- xts(df.pi[,-1], order.by = df.pi$date)

pi.xts <- log(pi.xts)

plot(unemployment.xts)

plot(realGDP.xts)

plot(nominalGDP.xts)

plot(GDPdeflator.xts)

plot(pi.xts)

Detrending the data as follows:

lin.modrGDP <- lm(realGDP.xts ~ time(realGDP.xts)) ## construct a linear model with the ts 
linear.rGDP <- xts(lin.modrGDP$fitted.values, order.by = df.pi$date) ## extract the trend line
cyclical.rGDP <- realGDP.xts-linear.rGDP ## subtract the trend line from the original ts object
plot(cyclical.rGDP) ##Boom, cyclical variations of real GDP, detrended.

ts.rGDP <- ts_ts(cyclical.rGDP)

plot(ts.rGDP)

rGDP.stl_decomp <- stl(ts.rGDP, s.window = "periodic")

rGDP.seas_comp <- rGDP.stl_decomp$time.series[, "seasonal"]

plot(rGDP.seas_comp)

    deseas.rGDP <- ts.rGDP - rGDP.seas_comp
    
    plot(deseas.rGDP)

adf.rGDP <- ur.df(deseas.rGDP, type = "none", selectlags = c("AIC")) ##Augmented Dickey-Fuller test for stationarity after we've removed the trend and the intercept
summary(adf.rGDP)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.138063 -0.002379  0.000746  0.004129  0.035587 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.31990    0.08351  -3.831 0.000207 ***
## z.diff.lag -0.30510    0.08765  -3.481 0.000703 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01555 on 117 degrees of freedom
## Multiple R-squared:  0.3024, Adjusted R-squared:  0.2905 
## F-statistic: 25.36 on 2 and 117 DF,  p-value: 7.077e-10
## 
## 
## Value of test-statistic is: -3.8308 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

What we’ve done here is run an Augmented Dickey-Fuller test to test for stationarity in the detrended, log transformed real GDP data. The value of the test statistic is -4.9934, which is less than the critical values at the 1%, 5% and 10% levels, which strongly implies real GDP is stationary.

##The concept remains the same, so we repeat the process for the rest of our NatAcc data
## Nominal GDP

lin.modnGDP <- lm(nominalGDP.xts ~ time(nominalGDP.xts))
linear.nGDP <- xts(lin.modnGDP$fitted.values, order.by = df.pi$date)
cyclical.nGDP <- nominalGDP.xts-linear.nGDP
plot(cyclical.nGDP)

ts.nGDP <- ts_ts(cyclical.nGDP)

plot(ts.nGDP)

nGDP.stl_decomp <- stl(ts.nGDP, s.window = "periodic")

nGDP.seas_comp <- nGDP.stl_decomp$time.series[, "seasonal"]

    deseas.nGDP <- ts.nGDP - nGDP.seas_comp
    
    plot(nGDP.seas_comp)

    plot(deseas.nGDP)

adf.nGDP <- ur.df(deseas.nGDP, type = "none", selectlags = c("AIC"))
summary(adf.nGDP)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.144689 -0.004462  0.001290  0.004667  0.044011 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.24043    0.07326  -3.282 0.001360 ** 
## z.diff.lag -0.31921    0.08749  -3.649 0.000395 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.017 on 117 degrees of freedom
## Multiple R-squared:  0.2609, Adjusted R-squared:  0.2482 
## F-statistic: 20.65 on 2 and 117 DF,  p-value: 2.093e-08
## 
## 
## Value of test-statistic is: -3.2818 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

The seasonaly adjusted nominal GDP series is stationary as well, with a test statistic on -4.7688, which is well below the 1pct. critical value, suggesting the datasetis stationary.

###################################################################################
## Unemloyment

lin.modUnem <- lm(unemployment.xts ~ time(unemployment.xts))
linear.Unem <- xts(lin.modUnem$fitted.values, order.by = df.pi$date)
cyclical.unem <- unemployment.xts-linear.Unem
plot(cyclical.unem)

ts.unem <- ts_ts(cyclical.unem)

plot(ts.unem)

unem.stl_decomp <- stl(ts.unem, s.window = "periodic")

unem.seas_comp <- unem.stl_decomp$time.series[, "seasonal"]

plot(unem.seas_comp)

    deseas.unem <- ts.unem - unem.seas_comp
    
    plot(deseas.unem)

adf.unem <- ur.df(deseas.unem, type = "none", selectlags = c("AIC"))
summary(adf.unem)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.196568 -0.010125 -0.001771  0.007795  0.096243 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)   
## z.lag.1    -0.16209    0.05774  -2.807  0.00586 **
## z.diff.lag -0.29106    0.08771  -3.318  0.00121 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02678 on 117 degrees of freedom
## Multiple R-squared:  0.1909, Adjusted R-squared:  0.1771 
## F-statistic: 13.81 on 2 and 117 DF,  p-value: 4.136e-06
## 
## 
## Value of test-statistic is: -2.8072 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

SImilar to the real GDP ADF test, our unemployment data is stationary since the test statistic of -4.7899 is less than (more negative than) the critical values at the 1% level (-2.6), the 5% level (-1.95) and the 10% level (-1.61). This means we can reject the null hypothesis and our unemployment data is stationary, which means it can be used in a VAR model as is.

###################################################################################

lin.modDeflator <- lm(GDPdeflator.xts ~ time(GDPdeflator.xts))
linear.Deflator <- xts(lin.modDeflator$fitted.values, order.by = df.pi$date)
cyclical.Deflator <- GDPdeflator.xts-linear.Deflator
plot(cyclical.Deflator)

ts.def <- ts_ts(cyclical.Deflator)

plot(ts.def)

def.stl_decomp <- stl(ts.def, s.window = "periodic")

def.seas_comp <- def.stl_decomp$time.series[, "seasonal"]

plot(def.seas_comp)

    deseas.def <- ts.def - def.seas_comp
    
    plot(deseas.def)

adf.def <- ur.df(deseas.def, type = "none", selectlags = c("AIC"))
summary(adf.def)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0095728 -0.0023548 -0.0002321  0.0017629  0.0118830 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.11791    0.05206  -2.265 0.025366 *  
## z.diff.lag -0.30853    0.08798  -3.507 0.000644 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003747 on 117 degrees of freedom
## Multiple R-squared:  0.1714, Adjusted R-squared:  0.1572 
## F-statistic:  12.1 on 2 and 117 DF,  p-value: 1.671e-05
## 
## 
## Value of test-statistic is: -2.2648 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

The Deflator data exhibits the same non-stationarity of the nominal GDP data. Which is implied, again, by the fact that the test statistic lies between the 5% and the 10% critical values, thus we fail to reject the null hypothesis here as well. Thus, the deflator is somewhat unsuitable for the analysis of monetary policy effects on inflation (as is) and I currently have no clue how to make it so, without manipulating the data too much and essentially “t” hacking).

Since we still require some sort of measure of inflation, we may be forced to consider using CPI data, which I was avoiding because I had already included CPI in calculating the user cost of money in the our törnqvist index. My intuition is that this “double counting” may lead to correlations which are somewhat spurious.

## We haven't detrended the inflation data lol.


lin.modPi <- lm(pi.xts ~ time(pi.xts))
linear.pi <- xts(lin.modPi$fitted.values, order.by = df.pi$date)
cyclical.Pi <- pi.xts-linear.pi
plot(cyclical.Pi)

ts.pi <- ts_ts(cyclical.Pi)

plot(ts.pi)

pi.stl_decomp <- stl(ts.pi, s.window = "periodic")

pi.seas_comp <- pi.stl_decomp$time.series[, "seasonal"]

plot(pi.seas_comp)

    deseas.pi <- ts.pi - pi.seas_comp
    
    plot(deseas.pi)

adf.inf <- ur.df(cyclical.Pi, type = "none", selectlags = c("AIC"))
summary(adf.inf)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52028 -0.04773  0.00802  0.07682  0.33237 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.15072    0.04249  -3.547 0.000561 ***
## z.diff.lag  0.33347    0.08776   3.800 0.000231 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1357 on 117 degrees of freedom
## Multiple R-squared:  0.1555, Adjusted R-squared:  0.141 
## F-statistic: 10.77 on 2 and 117 DF,  p-value: 5.091e-05
## 
## 
## Value of test-statistic is: -3.5472 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

After detrending the log transformed CPI data, we realise that this measure of inflation is more likely to be stationary than the deflator measure of prices. Though we fail to reject the null hypothesis (non-stationarity) at the 1% level, we can reject it at the 5% level, which means our wrangled CPI time series can be treated as being stationary. ——————————————————————————– ## Money Supply construction

The next step is to construct our index number for the money supply.

This part is a bit more involved. We begin by deciding what our benchmark rate is going to be.

For brevity’s sake we’ll set it as the Bank-Rate.

df.bankrate <- df.MonthlyQB %>% ## our benchmark rate in the formula for user cost: [(R-irr)/(1+R)]
  dplyr::select(
    date,
    KBP1401M ## Bankrate
  )

df.Yield <- df.QBweekly%>%
  dplyr::select(
    date,
    KBP1466W
  )

yield <- xts(df.Yield[,-1], order.by = df.QBweekly$date)

mo.yield <- to.monthly(yield, indexAt = "2013-01-01", drop.time = TRUE)



bankrate.xts <- xts(df.bankrate[,-1], order.by = df.bankrate$date)%>%
  na_ma()

plot(bankrate.xts)

align.benchmark <- bankrate.xts[index(irr.xts)]## aligning the indexes of the benchmark rate data

From here we take each observation of the benchmark rate and turn it into a vector of values. This vector is then used to calculate the user cost of money for each asset included in our monetary aggregate. The user cost is = (rb-or)/rb+1, where rb is the benchmark rate, and or is the own rate of return for the asset, which is logged in the irr time series object.

Using a for loop, we can calculate the user cost across the irr columns (which represent asset rates).

benchmrk.vec <- as.numeric(align.benchmark) ## create a numeric vector of the indexed bank rate data

usercost.xts <- irr.xts
for(i in 1:ncol(irr.xts)) { ## cycle through each column of the xts object
  asset_rate <- as.numeric(irr.xts[, i]) ## treating each irr entry as a numeric object
  usercost.xts[, i] <- (benchmrk.vec - asset_rate) / (1 + benchmrk.vec) ## calculate the user cost using irr and benchmark rate 
}

plot(usercost.xts)

The expenditure shares of each asset in each period can be calculated as s_{i,t} = (user cost{i,t} * quantity{i,t}) / total_expenditure_t. Here it becomes clear that the user cost plays the roll of the price of holding a particular asset as money. The a row-wise sum of expenditure shares should add up to one in every time period.

common_dates <- index(MonQuant.xts) ##alining the indexes
usercost.xts <- usercost.xts[common_dates,]
quant.money.xts <- MonQuant.xts[common_dates, ]
 


#  Calculate expenditure shares (s_{i,t})
#    s_{i,t} = (user_cost_{i,t} * quantity_{i,t}) / total_expenditure_t


expenditure.xts <- usercost.xts * quant.money.xts
total_expenditure <- rowSums(expenditure.xts, na.rm = TRUE)
shares <- expenditure.xts / total_expenditure ## s_{i,t}


plot(shares)

plot(expenditure.xts)

These are the expendature shares over time. Yes, the expendature is negative, that,s because it is user cost is negative

# 3. Compute lagged quantities and shares
lag_quantities <- lag.xts(quant.money.xts, 1)
lag_shares <- lag.xts(shares, 1)

# 4. Calculate Törnqvist index components
#    a. Quantity growth ratios
lag_quantities[lag_quantities < 1e-10] <- NA
growth_ratios <- quant.money.xts / lag_quantities
growth_ratios[growth_ratios <= 0 | is.infinite(growth_ratios)] <- NA
growth_ratios[growth_ratios > 1e100] <- NA

#    b. Weighted average of shares (t and t-1), this is the chaining of weights
weights <- (shares + lag_shares) / 2

#    c. Compute log-growth contributions
log_growth <- log(growth_ratios) * weights

# 5. Sum contributions across assets (for each period)
sum_log_growth <- rowSums(log_growth, na.rm = TRUE)

sumlog.growth.xts <- xts(sum_log_growth, order.by = common_dates)

summary(sumlog.growth.xts)  # Look for -Inf/Inf

##      Index            sumlog.growth.xts  
##  Min.   :2013-01-01   Min.   :-0.033261  
##  1st Qu.:2015-07-01   1st Qu.:-0.002606  
##  Median :2018-01-01   Median : 0.004283  
##  Mean   :2017-12-31   Mean   : 0.004690  
##  3rd Qu.:2020-07-01   3rd Qu.: 0.011814  
##  Max.   :2023-01-01   Max.   : 0.034683

cumulative_log_growth <- cumsum(replace(sum_log_growth, is.na(sum_log_growth), 0))
tornqvist_vals <- 100 * exp(cumulative_log_growth - cumulative_log_growth[1])
tornqvist_index <- xts(tornqvist_vals, order.by = common_dates)

n <- length(common_dates)
tornqvist_vals <- rep(NA_real_, n)
tornqvist_vals[1] <- 100  # Base period

for (i in 2:n) {
  if (!is.na(sum_log_growth[i])) {
    tornqvist_vals[i] <- tornqvist_vals[i-1] * exp(sum_log_growth[i])
  }
}
# 8. Create final xts object
tornqvist_index <- xts(tornqvist_vals, order.by = common_dates)



# 7. (Optional) Rename and inspect
colnames(tornqvist_index) <- "Tornqvist_Index"
head(tornqvist_index)

##            Tornqvist_Index
## 2013-01-01       100.00000
## 2013-02-01        99.83094
## 2013-03-01        99.86428
## 2013-04-01        99.81844
## 2013-05-01       101.76989
## 2013-06-01        99.99276

plot(sumlog.growth.xts)

plot(tornqvist_index)

Mt <- ts_ts(tornqvist_index)

plot(Mt)

Mt.stl_decomp <- stl(Mt, s.window = "periodic")

Mt.seas_comp <- Mt.stl_decomp$time.series[, "seasonal"]

    deseas.Mt <- Mt - Mt.seas_comp
    
    plot(deseas.Mt)

We then de-trend the log transformation of the data using a linear model as we did with the National Accounts Data:

logMt <- log(Mt)

lin.modMt <- lm(logMt~ time(Mt))
linear.Mt <- xts(lin.modMt$fitted.values, order.by = df.pi$date)
lin.Mt <- ts_ts(linear.Mt)
cyclical.Mt <- logMt-lin.Mt
plot(cyclical.Mt)

Mt <- ts_ts(cyclical.Mt)

plot(Mt)

Mt.stl_decomp <- stl(Mt, s.window = "periodic")

Mt.seas_comp <- Mt.stl_decomp$time.series[, "seasonal"]

plot(Mt.seas_comp )

    deseas.Mt <- Mt - Mt.seas_comp
    
    plot(deseas.Mt)

adf.Mt <- ur.df(cyclical.Mt, type = "none", selectlags = c("AIC"))
summary(adf.Mt)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.031518 -0.007107 -0.000054  0.008172  0.031708 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)  
## z.lag.1    -0.09068    0.04355  -2.082   0.0395 *
## z.diff.lag -0.09389    0.09320  -1.007   0.3158  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01296 on 117 degrees of freedom
## Multiple R-squared:  0.05458,    Adjusted R-squared:  0.03842 
## F-statistic: 3.377 on 2 and 117 DF,  p-value: 0.03749
## 
## 
## Value of test-statistic is: -2.0825 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

The detrended Mt, is stationary, as the test statistic is -2.3233, which is less than the critical value at 5%, which allows us to reject the null hypothesis.

From here we have the tough job of running VARs on this data to come up with a characteristic equation that can be turned into Structural Auto-regressive models.

The first step is to bind our relevant time series together into an xts object with three variables. I chose three because it is only real GDP, CPI, Unemployment and Mt that are likely to be stationary according to the ADF test.

Herein comes the question; to what extent does monetary policy influence factors like real GDP, CPI and the rate of unemployment.

So this means our data will need to be ordered in a specific manner in order for the structural Autoregression to make economic sense.

I believe the following ordering is most in line with economic theory (rGDP, Unem, Mt, CPI). A Cholesky decomposition on this ordering implies:

Only changes in real output can affect real output contemporaneously.
Changes on real output and unemployment can affect unemployment contemporaneously.
Changes in real GDP, Unemployment and the money supply can affect the money supply contemporaneously.
All variables can affect the rate of inflation.

Point number 3 may require some justification. Though I believe the explanation to be straightforward.

A shock to real income is likely to immediately influence the demand to hold money in different forms (perhaps a movement of money from less liquid accounts such as savings accounts to more liquid accounts such as check accounts).

A shock to Unemployment, similarly influences the effective demand to hold money as these represent changes to the level of an agents permanent income which, according to the permanent income hypothesis, alters the optimal the saving and spending behaviour of agents. Since it is the through money that spending and saving takes place, then the demand for holding money will react to these aforementioned shocks.

point 4 is somewhat of a sticking point as well but it is also justifiable with economic theory. This can be shown by constructing an Aggregate Supply and Demand schedule, and noting that the output axis is obviously contemporaneous with inflation as that is on the y axis. The negative output gap, which is proportional to the rate of unemployment is also reflected in the same diagram and on the same footing as changes in real output (changes in real output imply changes in resource employment). Though the output gap represents spare capacity of all kinds, not just unused labour, we assume that, in the short term, it is through changes in labour employment that represents the most income elastic resource for firms and businesses. ———————————————————————————————

cycYUMP <- cbind(deseas.rGDP, deseas.unem, deseas.Mt, deseas.pi)## stationary version
colnames(cycYUMP) <- c("rGDP-dt", "Unem-dt", "Mt-dt", "pi-dt")
plot(cycYUMP)

The YUMP dataset is just the Log transformed Data, while cycYUMP represents the detrended data.

info.var <- vars::VARselect(cycYUMP, lag.max = 12, type = "none")
cycYump.Lag <- info.var$selection["AIC(n)"]

cycYUMP.model <- vars::VAR(
  y = cycYUMP,
  p = cycYump.Lag,
  type = "none",
  season = NULL,
  exogen = NULL
)

summary(cycYUMP.model)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: rGDP.dt, Unem.dt, Mt.dt, pi.dt 
## Deterministic variables: none 
## Sample size: 118 
## Log Likelihood: 1162.597 
## Roots of the characteristic polynomial:
## 0.9168 0.8735 0.8735 0.7814 0.5755 0.5755 0.5299 0.5299 0.4278 0.4278 0.1678 0.1016
## Call:
## vars::VAR(y = cycYUMP, p = cycYump.Lag, type = "none", exogen = NULL)
## 
## 
## Estimation results for equation rGDP.dt: 
## ======================================== 
## rGDP.dt = rGDP.dt.l1 + Unem.dt.l1 + Mt.dt.l1 + pi.dt.l1 + rGDP.dt.l2 + Unem.dt.l2 + Mt.dt.l2 + pi.dt.l2 + rGDP.dt.l3 + Unem.dt.l3 + Mt.dt.l3 + pi.dt.l3 
## 
##             Estimate Std. Error t value Pr(>|t|)    
## rGDP.dt.l1  0.198498   0.147532   1.345   0.1813    
## Unem.dt.l1 -0.046585   0.092153  -0.506   0.6142    
## Mt.dt.l1    0.220563   0.138978   1.587   0.1155    
## pi.dt.l1    0.050792   0.009826   5.169 1.11e-06 ***
## rGDP.dt.l2  0.264864   0.164136   1.614   0.1096    
## Unem.dt.l2  0.050999   0.104902   0.486   0.6279    
## Mt.dt.l2   -0.053987   0.167374  -0.323   0.7477    
## pi.dt.l2   -0.025063   0.015846  -1.582   0.1167    
## rGDP.dt.l3  0.226986   0.157340   1.443   0.1521    
## Unem.dt.l3 -0.007163   0.090913  -0.079   0.9374    
## Mt.dt.l3   -0.126213   0.135698  -0.930   0.3544    
## pi.dt.l3   -0.019818   0.011419  -1.735   0.0856 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.0137 on 106 degrees of freedom
## Multiple R-Squared: 0.5448,  Adjusted R-squared: 0.4933 
## F-statistic: 10.57 on 12 and 106 DF,  p-value: 1.823e-13 
## 
## 
## Estimation results for equation Unem.dt: 
## ======================================== 
## Unem.dt = rGDP.dt.l1 + Unem.dt.l1 + Mt.dt.l1 + pi.dt.l1 + rGDP.dt.l2 + Unem.dt.l2 + Mt.dt.l2 + pi.dt.l2 + rGDP.dt.l3 + Unem.dt.l3 + Mt.dt.l3 + pi.dt.l3 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## rGDP.dt.l1 -0.44885    0.23810  -1.885  0.06215 .  
## Unem.dt.l1  0.46154    0.14872   3.103  0.00245 ** 
## Mt.dt.l1    0.40328    0.22430   1.798  0.07502 .  
## pi.dt.l1    0.07805    0.01586   4.922 3.16e-06 ***
## rGDP.dt.l2 -0.03583    0.26490  -0.135  0.89266    
## Unem.dt.l2  0.31887    0.16930   1.883  0.06238 .  
## Mt.dt.l2   -0.15432    0.27012  -0.571  0.56900    
## pi.dt.l2   -0.04484    0.02557  -1.753  0.08245 .  
## rGDP.dt.l3 -0.21799    0.25393  -0.858  0.39257    
## Unem.dt.l3  0.05736    0.14672   0.391  0.69663    
## Mt.dt.l3   -0.22038    0.21900  -1.006  0.31657    
## pi.dt.l3   -0.02430    0.01843  -1.318  0.19025    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.02211 on 106 degrees of freedom
## Multiple R-Squared: 0.7793,  Adjusted R-squared: 0.7543 
## F-statistic: 31.18 on 12 and 106 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation Mt.dt: 
## ====================================== 
## Mt.dt = rGDP.dt.l1 + Unem.dt.l1 + Mt.dt.l1 + pi.dt.l1 + rGDP.dt.l2 + Unem.dt.l2 + Mt.dt.l2 + pi.dt.l2 + rGDP.dt.l3 + Unem.dt.l3 + Mt.dt.l3 + pi.dt.l3 
## 
##             Estimate Std. Error t value Pr(>|t|)    
## rGDP.dt.l1  0.090472   0.101859   0.888    0.376    
## Unem.dt.l1 -0.084820   0.063625  -1.333    0.185    
## Mt.dt.l1    0.659297   0.095954   6.871 4.53e-10 ***
## pi.dt.l1    0.008244   0.006784   1.215    0.227    
## rGDP.dt.l2  0.113783   0.113323   1.004    0.318    
## Unem.dt.l2 -0.032440   0.072427  -0.448    0.655    
## Mt.dt.l2    0.166487   0.115559   1.441    0.153    
## pi.dt.l2   -0.004109   0.010940  -0.376    0.708    
## rGDP.dt.l3  0.023282   0.108631   0.214    0.831    
## Unem.dt.l3  0.039172   0.062769   0.624    0.534    
## Mt.dt.l3   -0.020880   0.093690  -0.223    0.824    
## pi.dt.l3    0.003598   0.007884   0.456    0.649    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.009458 on 106 degrees of freedom
## Multiple R-Squared: 0.8937,  Adjusted R-squared: 0.8816 
## F-statistic: 74.25 on 12 and 106 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation pi.dt: 
## ====================================== 
## pi.dt = rGDP.dt.l1 + Unem.dt.l1 + Mt.dt.l1 + pi.dt.l1 + rGDP.dt.l2 + Unem.dt.l2 + Mt.dt.l2 + pi.dt.l2 + rGDP.dt.l3 + Unem.dt.l3 + Mt.dt.l3 + pi.dt.l3 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## rGDP.dt.l1  0.35978    1.44018   0.250   0.8032    
## Unem.dt.l1 -1.52379    0.89958  -1.694   0.0932 .  
## Mt.dt.l1    1.29218    1.35668   0.952   0.3430    
## pi.dt.l1    1.20415    0.09592  12.554   <2e-16 ***
## rGDP.dt.l2 -0.53768    1.60227  -0.336   0.7379    
## Unem.dt.l2  0.61302    1.02404   0.599   0.5507    
## Mt.dt.l2   -0.57560    1.63389  -0.352   0.7253    
## pi.dt.l2   -0.29132    0.15469  -1.883   0.0624 .  
## rGDP.dt.l3 -1.18303    1.53593  -0.770   0.4429    
## Unem.dt.l3  0.83042    0.88749   0.936   0.3516    
## Mt.dt.l3   -0.86081    1.32467  -0.650   0.5172    
## pi.dt.l3   -0.04470    0.11147  -0.401   0.6893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.1337 on 106 degrees of freedom
## Multiple R-Squared: 0.8254,  Adjusted R-squared: 0.8057 
## F-statistic: 41.77 on 12 and 106 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           rGDP.dt   Unem.dt      Mt.dt      pi.dt
## rGDP.dt 1.876e-04 2.454e-04  8.324e-07  8.771e-05
## Unem.dt 2.454e-04 4.877e-04  9.579e-06  6.242e-05
## Mt.dt   8.324e-07 9.579e-06  8.928e-05 -9.666e-05
## pi.dt   8.771e-05 6.242e-05 -9.666e-05  1.788e-02
## 
## Correlation matrix of residuals:
##          rGDP.dt Unem.dt     Mt.dt    pi.dt
## rGDP.dt 1.000000 0.81124  0.006431  0.04788
## Unem.dt 0.811241 1.00000  0.045907  0.02114
## Mt.dt   0.006431 0.04591  1.000000 -0.07650
## pi.dt   0.047882 0.02114 -0.076499  1.00000

Hopefully it is clear that the de-trended data and the simple log transformed data are somewhat similar in their implications, particularly with regard to changes in inflation where some lagged values of the money supply, real GDP and inflation itself are significant at a level smaller than 10%.

From here we move towards a structural VAR. We will begin by constructing our contemporaneous effects matrix.

a.mat<- diag(4)
diag(a.mat)<- NA
a.mat[2, 1] <- NA
a.mat[3, 1] <- NA
a.mat[3, 2] <- NA
a.mat[4, 1] <- NA
a.mat[4, 2] <- NA
a.mat[4, 3] <- NA

print(a.mat)

##      [,1] [,2] [,3] [,4]
## [1,]   NA    0    0    0
## [2,]   NA   NA    0    0
## [3,]   NA   NA   NA    0
## [4,]   NA   NA   NA   NA

We have imposed 6 restrictions on the upper triangle which satisfies the K(K-1)/2 when K=4.

For our B matrix (For individual shocks)

b.mat <- diag(4)
diag(b.mat)<- NA
print(b.mat)

##      [,1] [,2] [,3] [,4]
## [1,]   NA    0    0    0
## [2,]    0   NA    0    0
## [3,]    0    0   NA    0
## [4,]    0    0    0   NA

svarcycYUMP <- vars::SVAR(cycYUMP.model, Amat = a.mat, Bmat = b.mat, max.iter = 1000, hessian = TRUE)

## Warning in vars::SVAR(cycYUMP.model, Amat = a.mat, Bmat = b.mat, max.iter =
## 1000, : The AB-model is just identified. No test possible.

svarcycYUMP

## 
## SVAR Estimation Results:
## ======================== 
## 
## 
## Estimated A matrix:
##          rGDP.dt  Unem.dt Mt.dt pi.dt
## rGDP.dt  1.00000  0.00000 0.000     0
## Unem.dt -1.30627  1.00000 0.000     0
## Mt.dt    0.05644 -0.04698 1.000     0
## pi.dt   -0.82209  0.26676 1.058     1
## 
## Estimated B matrix:
##         rGDP.dt Unem.dt    Mt.dt  pi.dt
## rGDP.dt  0.0137 0.00000 0.000000 0.0000
## Unem.dt  0.0000 0.01298 0.000000 0.0000
## Mt.dt    0.0000 0.00000 0.009438 0.0000
## pi.dt    0.0000 0.00000 0.000000 0.1331

Our Impulse Response functions will follow from here.

YU_results <- vars::irf(cycYUMP.model, 
                   impulse = "rGDP.dt", 
                   response = "Unem.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(YU_results)

YM_results <- vars::irf(cycYUMP.model, 
                   impulse = "rGDP.dt", 
                   response = "Mt.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(YM_results)

MU_results <- vars::irf(cycYUMP.model, 
                   impulse = "Mt.dt", 
                   response = "Unem.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(MU_results)

MP_results <- vars::irf(cycYUMP.model, 
                   impulse = "Mt.dt", 
                   response = "pi.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(MP_results)

YP_results <- vars::irf(cycYUMP.model, 
                   impulse = "rGDP.dt", 
                   response = "pi.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(YP_results)

UP_results <- vars::irf(cycYUMP.model, 
                   impulse = "Unem.dt", 
                   response = "pi.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(UP_results)

PU_results <- vars::irf(cycYUMP.model, 
                   impulse = "pi.dt", 
                   response = "Unem.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(PU_results)

UY_results <- vars::irf(cycYUMP.model, 
                   impulse = "Unem.dt", 
                   response = "rGDP.dt",
                   n.ahead = 40, 
                   ortho = TRUE, 
                   boot = TRUE, 
                   ci = 0.95)
plot(UY_results)

## For longer term relationships

cyc.lt.YUMP <- vars::BQ(cycYUMP.model)

summary(cyc.lt.YUMP)

## 
## SVAR Estimation Results:
## ======================== 
## 
## Call:
## vars::BQ(x = cycYUMP.model)
## 
## Type: Blanchard-Quah 
## Sample size: 118 
## Log Likelihood: 1137.288 
## 
## Estimated contemporaneous impact matrix:
##           rGDP.dt    Unem.dt     Mt.dt     pi.dt
## rGDP.dt 0.0120628  0.0014002 -0.004497 -0.004469
## Unem.dt 0.0142397  0.0146662 -0.005051 -0.006740
## Mt.dt   0.0006908 -0.0001864  0.007400 -0.005847
## pi.dt   0.0657330  0.0067860  0.059798  0.099705
## 
## Estimated identified long run impact matrix:
##          rGDP.dt  Unem.dt   Mt.dt  pi.dt
## rGDP.dt  0.05294  0.00000 0.00000 0.0000
## Unem.dt -0.12689  0.08638 0.00000 0.0000
## Mt.dt    0.11237 -0.03410 0.05358 0.0000
## pi.dt   -0.09351  0.03612 0.39487 0.7561
## 
## Covariance matrix of reduced form residuals (*100):
##           rGDP.dt   Unem.dt      Mt.dt     pi.dt
## rGDP.dt 1.877e-02 0.0245145  9.266e-05  0.008790
## Unem.dt 2.451e-02 0.0488811  9.131e-04  0.006148
## Mt.dt   9.266e-05 0.0009131  8.946e-03 -0.009629
## pi.dt   8.790e-03 0.0061476 -9.629e-03  1.788371

irf.rGDP <- irf(cyc.lt.YUMP, impulse = "rGDP.dt", boot = FALSE, n.ahead = 40)
irf.Unem <- irf(cyc.lt.YUMP, impulse = "Unem.dt", boot = FALSE, n.ahead = 40)
irf.Mt <- irf(cyc.lt.YUMP, impulse = "Mt.dt", boot = FALSE, n.ahead = 40)
irf.Pi <- irf(cyc.lt.YUMP, impulse = "pi.dt", boot = FALSE, n.ahead = 40)

plot(irf.rGDP)

plot(irf.Unem)

plot(irf.Mt)

plot(irf.Pi)

Time Series Analysis of Monetary Policy

K.D. Leshaba

2025-09-19

Goals

Stationarity makes progress

Samadb and retrieving data

Pause, real variables.