Expectation-augmented Phillips curve.Rmd

---
title: "Expectation-augmented Phillips curve"
author: "Kyle Deng"
date: "Last updated at `r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    df_print: paged
---

```{r Knitr_Global_Options, include=FALSE, cache=FALSE}
library(knitr)
library(tidyverse)
library(lubridate)
library(dotenv)
library(fredr)
library(ggrepel)
library(ggfortify)
library(xts)
library(AER)
library(arsenal)
library(stargazer)
opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, 
               autodep = TRUE, tidy = FALSE, cache = TRUE,
               fig.dim=c(7.7,4.5))
load_dot_env("cred.env")
# Cpi series
cpiaucsl <- fredr("CPIAUCSL", frequency = "q") %>%
  na.omit() %>%
  select(date, value) %>%
  rename(CPIAUCSL = value)
# Core cpi series
cpilfesl <- fredr("CPILFESL", frequency = "q") %>%
  na.omit() %>%
  select(date, value) %>%
  rename(CPILFESL = value)
# Unemployment rate
unrate <- fredr("UNRATE", frequency = "q") %>%
  na.omit() %>%
  select(date, value) %>%
  rename(UNRATE = value)
# Natural rate of unemployment 
nroust <- fredr("NROUST", frequency = "q", observation_end = as.Date("2021-04-01")) %>%
  na.omit() %>%
  select(date, value) %>%
  rename(NROUST = value)
# Joining the data set together
data <- left_join(cpiaucsl,cpilfesl, by = "date") %>%
  left_join(unrate, by = "date") %>%
  left_join(nroust, by = "date")
```


#### Obtain the inflation surprises and unemployment gap according to Coiban and Gorodnichenko(2015).

##### Backward-looking inflation expecatation formula:
$$E_t\pi_{t+1}=\frac{1}{4}(\pi_{t-1}+\pi_{t-2}+\pi_{t-3}+\pi_{t-4})$$

##### Inflation surprises:
$$\pi_t-E_t\pi_{t+1}$$

##### Unemployment gap:
$$UE_t^{gap}=UE_t-UE_t^{natural}$$
```{r}
data$cpiaucsl_a = (data$CPIAUCSL/lag(data$CPIAUCSL, 4) - 1) * 100 
data$cpiaucsl_q = (data$CPIAUCSL/lag(data$CPIAUCSL, 1) - 1) * 400
data$back_cpiaucsl_qa = data$cpiaucsl_q - lag(data$cpiaucsl_a, 1)
data$cpilfesl_a = (data$CPILFESL/lag(data$CPILFESL, 4) - 1) * 100 
data$cpilfesl_q = (data$CPILFESL/lag(data$CPILFESL, 1) - 1) * 400
data$back_cpilfesl_qa = data$cpilfesl_q - lag(data$cpilfesl_a, 1)
data$un_gap <- data$UNRATE - data$NROUST

# Delete observation before 1985Q1  
# create a split point base on the cut-off point chosen in Figure 5 of Coiban and Gorodnichenko (2015) for my initial analysis
data <- data[-c(1:which(data$date == "1984-10-01")),] %>%
  mutate(label1 = case_when(
    date < "2007-10-01"  ~ "1985Q1-2007Q3",
    date >= "2007-10-01" ~ "2007Q4-2021Q2"
  ))
# Create a time series stamp
data$t <- seq.int(nrow(data))
```

### Summary statistics
```{r}
table1::label(data$back_cpiaucsl_qa) <- "Inflation surprises"
table1::label(data$back_cpilfesl_qa) <- "Core inflation surprises"
table1::label(data$un_gap) <- "Unemployment gap"
table1::table1(~back_cpiaucsl_qa + back_cpilfesl_qa + un_gap | label1, data = data)
```


### Time series plots
```{r}
date <- seq(data$date[1], data$date[146], by = '+3 month')
plot(xts(data$back_cpiaucsl_qa, order.by = date), main = "CPI inflation surprise")
plot(xts(data$back_cpilfesl_qa, order.by = date), main = "Core CPI inflation surprise")
plot(xts(data$un_gap, order.by = date), main = "Unemployment gap")
```

### Scatter plot to show the Phillips curve relationship
#### Phillips curve with regular CPI inflation
```{r}
ggplot() +
  geom_point(data = data, mapping = aes(x = un_gap, y = back_cpiaucsl_qa, color = label1, shape = label1)) + 
  geom_smooth(data = data, mapping = aes(x = un_gap, y = back_cpiaucsl_qa), method = "lm", color = "black", se = F, size = 0.75) +
  coord_cartesian(ylim = c(-5.5, 5), xlim = c(-1.5, 5)) +
  xlab("Unemployment gap") +
  ylab("CPI Inflation surprises") 
```

2008Q4 and 2020Q2 are out of bound so are not reported in the regular CPI inflation graph

#### Phillips curve with core CPI inflation
```{r}
ggplot() +
  geom_point(data = data, mapping = aes(x = un_gap, y = back_cpilfesl_qa, color = label1, shape = label1)) + 
  geom_smooth(data = data, mapping = aes(x = un_gap, y = back_cpilfesl_qa), method = "lm", color = "black", se = F, size = 0.75) +
  coord_cartesian(ylim = c(-4, 3), xlim = c(-1.5, 9)) +
  xlab("Unemployment gap") +
  ylab("Core CPI Inflation surprises") 
```

### OLS and IV estimate of the Phillips curve relationship for the full sample
##### OLS formula:
$$\pi_t-E_t\pi_{t+1}=c+\kappa UE_t^{gap}+v_t$$

#### Instrumental variable, two-stage least squares formula:
##### First stage regression:
$$\widehat{UE_t^{gap}}=\alpha+\beta UE_{t-1}^{gap}+\epsilon_t$$

##### Outcome regression:
$$\pi_t-E_t\pi_{t+1}=c+\kappa \widehat{UE_{t-1}^{gap}}+u_t$$
```{r}
## OLS Phillips curve model for regular CPI 
full_ols <- lm(back_cpiaucsl_qa ~ un_gap, data = data)
## OLS Phillips curve model for core CPI 
full_ols_core <- lm(back_cpilfesl_qa ~ un_gap, data = data)
## IV Phillips curve model for regular CPI 
full_iv <- ivreg(back_cpiaucsl_qa ~ un_gap | lag(un_gap, 1), data = data)
## IV Phillips curve model for core CPI 
full_iv_core <- ivreg(back_cpilfesl_qa ~ un_gap | lag(un_gap, 1), data = data)

# Newey-West standard errors
robust3 <- sqrt(diag(NeweyWest(full_ols, prewhite = F, adjust = T)))
robust4 <- sqrt(diag(NeweyWest(full_iv, prewhite = F, adjust = T)))
robust3c <- sqrt(diag(NeweyWest(full_ols_core, prewhite = F, adjust = T)))
robust4c <- sqrt(diag(NeweyWest(full_iv_core, prewhite = F, adjust = T)))
```

#### Estimates for the models with Newey-West standard errors
##### OLS Phillips curve model for regular CPI
```{r}
coeftest(full_ols, vcov.=NeweyWest(full_ols, prewhite = F, adjust = T))
```

##### IV Phillips curve model for regular CPI 
```{r}
coeftest(full_iv, vcov.=NeweyWest(full_iv, prewhite = F, adjust = T))
```

##### OLS Phillips curve model for core CPI 
```{r}
coeftest(full_ols_core, vcov.=NeweyWest(full_ols_core, prewhite = F, adjust = T))
```

##### IV Phillips curve model for core CPI 
```{r}
coeftest(full_iv_core, vcov.=NeweyWest(full_iv_core, prewhite = F, adjust = T))
```


### Exhaustive search for the split point
##### Chow-test statistics:
$$C=\frac{RSS-RSS_1-RSS_2}{RSS_1+RSS_2}\times\frac{T-2k}{k}$$
$k$ = number of parameters in the model

#### Analysis on the regular CPI inflation data
```{r}
RSS = sum(resid(full_iv)^2)

t <- ceiling(0.15 * nrow(data))
tt <- floor(0.85 * nrow(data))
chow <- 0
index <- 0
for(i in t:tt){
  data1 <- data %>% filter(t <= i)
  data2 <- data %>% filter(t > i)
  sample1 <- ivreg(back_cpiaucsl_qa ~ un_gap | lag(un_gap, 1) , data = data1)
  sample2 <- ivreg(back_cpiaucsl_qa ~ un_gap | lag(un_gap, 1) , data = data2)
  RSS1 = sum(resid(sample1)^2)
  RSS2 = sum(resid(sample2)^2)
  C = ((RSS - RSS1 - RSS2) / (RSS1 + RSS2)) * ((146 - 2*2) / 2)
  if(C > chow){
    chow = C
    index = i
  }
}
qalpha = qf(0.95,2,146-4)
ifelse(chow>qalpha, 
       "Decision: Reject Null of Stability with signif. level 5%",
       "Decision: Do NOT reject Null of Stability with signif. level 5%")
## Split point from the search
data$date[index]
## Create a dummy base on the split 
data$dummy = ifelse(data$date >= data$date[index], 1, 0)
## Create another label base on the split
data <- data %>% 
  mutate(label2 = case_when(
    date >= "1985-01-01" & date < "2008-07-01"   ~ "1985Q1-2008Q2",
    date >= "2008-07-01" ~ "2008Q3-2021Q2"
  ))
```

#### Analysis on the core CPI inflation data
```{r}
RSS_core = sum(resid(full_iv_core)^2)
chowc <- 0
indexc <- 0
for(i in t:tt){
  data1 <- data %>% filter(t <= i)
  data2 <- data %>% filter(t > i)
  sample1 <- ivreg(back_cpilfesl_qa ~ un_gap | lag(un_gap, 1) , data = data1)
  sample2 <- ivreg(back_cpilfesl_qa ~ un_gap | lag(un_gap, 1) , data = data2)
  RSS1c = sum(resid(sample1)^2)
  RSS2c = sum(resid(sample2)^2)
  Cc = ((RSS_core - RSS1c - RSS2c) / (RSS1c + RSS2c)) * ((146 - 2*2) / 2)
  if(Cc > chowc){
    chowc = Cc
    indexc = i
  }
}
ifelse(chowc>qalpha, 
       "Decision: Reject Null of Stability with signif. level 5%",
       "Decision: Do NOT reject Null of Stability with signif. level 5%")
## Split point from the search
data$date[indexc]
## Create a dummy base on the split 
data$dummyc = ifelse(data$date >= data$date[indexc], 1, 0)
## Create another label base on the split
data <- data %>% 
  mutate(label3 = case_when(
    date >= "1985-01-01" & date < "1991-01-01"   ~ "1985Q1-1990Q4",
    date >= "1991-01-01" ~ "1991Q1-2021Q2"
  ))
```
Although the Chow test indicates the core CPI inflation series is stable over the entire sample period, I still proceed with the split point to investigate the relationship out of curiosity.  

### Phillips curve on the regular CPI inflation data with the new split 
```{r}
data_pre = subset(data, label2 == "1985Q1-2008Q2")
data_post = subset(data, label2 == "2008Q3-2021Q2")

ggplot() +
  geom_point(data = data, mapping = aes(x = un_gap, y = back_cpiaucsl_qa, 
                                        color = label2, shape = label2)) +
  geom_text(data = data[142:146, ], aes(x = un_gap, y = back_cpiaucsl_qa, 
                                        label = date), position=position_dodge(width=1.5), 
                                    hjust = 1.5, size = 2.5) +
  geom_text_repel(data = data[142:146, ], aes(x = un_gap, y = back_cpiaucsl_qa, 
                                              label = date), size = 2.5) +
  geom_smooth(data = data_pre, mapping = aes(x = un_gap, y = back_cpiaucsl_qa), 
              method = "lm", se = F, color = "#F8766D", size = 1) +
  geom_smooth(data = data_post, mapping = aes(x = un_gap, y = back_cpiaucsl_qa), 
              method = 'lm', se = F, color = "#00BFC4", size = 1) +
  ggtitle("Time variation in the slope of the Phillips curve") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(ylim = c(-14, 5.5), xlim = c(-2, 8.5)) +
  xlab("Unemployment gap") +
  ylab("CPI inflation surprises") 
```

### Phillips curve on the core CPI inflation data with the new split
```{r}
data_prec = subset(data, label3 == "1985Q1-1990Q4")
data_postc = subset(data, label3 == "1991Q1-2021Q2")

ggplot() +
  geom_point(data = data, mapping = aes(x = un_gap, y = back_cpilfesl_qa, 
                                        color = label3, shape = label3)) +
  geom_text(data = data[142:146, ], aes(x = un_gap, y = back_cpilfesl_qa, label = date), 
            position=position_dodge(width=1.5), hjust = 1.5, size = 2.5) +
  geom_text_repel(data = data[142:146, ], aes(x = un_gap, y = back_cpilfesl_qa, 
                                              label = date), size = 2.5) +
  geom_smooth(data = data_prec, mapping = aes(x = un_gap, y = back_cpilfesl_qa), 
              method = "lm", se = F, color = "#F8766D", size = 1) +
  geom_smooth(data = data_postc, mapping = aes(x = un_gap, y = back_cpilfesl_qa), 
              method = 'lm', se = F, color = "#00BFC4", size = 1) +
  coord_cartesian(ylim = c(-4, 6.25), xlim = c(-1.75, 8.5)) +
  ggtitle("Time variation in the slope of the Phillips curve (Core CPI inflation)") +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Unemployment gap") +
  ylab("Core CPI inflation surprises") 
```

### Seperate IV regression on the 2 sub-periods to find out how much the slope has changed
````{r}
# Regular CPI inflation
split1 <- ivreg(back_cpiaucsl_qa ~ un_gap | un_gap + lag(un_gap, 1), data = data[1:index-1,])
# summary(split1)
robust1 <- sqrt(diag(NeweyWest(split1, prewhite = F, adjust = T)))
split2 <- ivreg(back_cpiaucsl_qa ~ un_gap | un_gap + lag(un_gap, 1), data = data[index:nrow(data),])
robust2 <- sqrt(diag(NeweyWest(split2, prewhite = F, adjust = T)))
# summary(split2)
```

##### 1985Q1-2008Q2 with rugular CPI inflation
```{r}
coeftest(split1, vcov.= NeweyWest(split1, prewhite = F, adjust = T))
```

##### 2008Q3-2021Q2 with rugular CPI inflation
```{r}
coeftest(split2, vcov.= NeweyWest(split2, prewhite = F, adjust = T))
```

```{r}
# Core CPI inflation
split1c <- ivreg(back_cpilfesl_qa ~ un_gap | lag(un_gap, 1), data = data[1:indexc-1,])
split2c <- ivreg(back_cpilfesl_qa ~ un_gap | lag(un_gap, 1), data = data[indexc:nrow(data),])
robust1c <- sqrt(diag(NeweyWest(split1c, prewhite = F, adjust = T)))
robust2c <- sqrt(diag(NeweyWest(split2c, prewhite = F, adjust = T)))
```

##### 1985Q1-1990Q4 with core CPI inflation
```{r}
coeftest(split1c, vcov.= NeweyWest(split1c, prewhite = F, adjust = T))
```

##### 1991Q1-2021Q2 with core CPI inflation
```{r}
coeftest(split2c, vcov.= NeweyWest(split2c, prewhite = F, adjust = T))
```


### Output the table of the regression results
##### Table for regular CPI inflation
```{r table1, results = "asis"}
stargazer(split1, split2, full_ols, full_iv, type = 'html', object.names = T, model.numbers = F, se = list(robust1, robust2, robust3, robust4), align = T, dep.var.labels = c("CPI inflation surprises"), covariate.labels = c("Unemployment gap hat"), header = F)
```


##### Table for core CPI inflation
```{r table2, results = "asis"}
stargazer(split1c, split2c, full_ols_core, full_iv_core, type = 'html', object.names = T, model.numbers = F, se = list(robust1c, robust2c, robust3c, robust4c), align = T, dep.var.labels = c("Core CPI inflation surprises"), covariate.labels = c("Unemployment gap hat"), header = F)
```