library(tidyverse)
library(forecast)
library(tseries)
library(zoo)
library(lubridate)
library(gridExtra)
library(knitr)
library(scales)

# Read raw without headers
raw <- read.csv("Prices-Inflation.csv", header = FALSE, stringsAsFactors = FALSE,
                na.strings = c("", "NA"))

# Row 1 = category names, Row 2 = sub-metric labels, Row 3+ = data
col_names <- c(
  "Date",
  "CPI_Index", "CPI_qq", "CPI_yy",
  "Tradable_Index", "Tradable_qq", "Tradable_yy",
  "NonTradable_Index", "NonTradable_qq", "NonTradable_yy",
  "Food_Index", "Food_qq", "Food_yy",
  "Housing_Index", "Housing_qq", "Housing_yy",
  "Transport_Index", "Transport_qq", "Transport_yy",
  "Petrol_Index", "Petrol_qq", "Petrol_yy",
  "CPI_ExPetrol_Index", "CPI_ExPetrol_qq", "CPI_ExPetrol_yy",
  "TrimmedMean_qq", "TrimmedMean_yy",
  "WeightedMedian_qq", "WeightedMedian_yy",
  "FactorModel_yy",
  "SectoralFactor_yy", "SectoralTradable_yy", "SectoralNonTradable_yy",
  "HPI_Index", "HPI_qq", "HPI_yy",
  "GDPDeflator_Index", "GDPDeflator_qq", "GDPDeflator_yy",
  "PPI_Input_Index", "PPI_Input_qq", "PPI_Input_yy",
  "PPI_Output_Index", "PPI_Output_qq", "PPI_Output_yy"
)

df <- raw[3:nrow(raw), 1:length(col_names)]
colnames(df) <- col_names

# Parse date and numeric columns
df <- df %>%
  mutate(Date = as.yearqtr(Date, format = "%b %Y")) %>%
  mutate(across(-Date, ~ suppressWarnings(as.numeric(.x)))) %>%
  arrange(Date) %>%
  filter(!is.na(CPI_Index))

cat("Data range:", format(min(df$Date)), "to", format(max(df$Date)), "\n")

cat("Rows:", nrow(df), " | Columns:", ncol(df), "\n")

kable(tail(df[, 1:7], 8), digits = 2, caption = "Recent CPI Data (last 8 quarters)")

# Use y/y% from 1989 (first full year of data)
cpi_yy <- df %>%
  filter(!is.na(CPI_yy)) %>%
  arrange(Date)

cpi_ts <- ts(cpi_yy$CPI_yy,
             start = c(year(min(as.Date(cpi_yy$Date))),
                       quarter(min(as.Date(cpi_yy$Date)))),
             frequency = 4)

cat("TS start:", start(cpi_ts), "| end:", end(cpi_ts),
    "| length:", length(cpi_ts), "\n")

stl_fit <- stl(cpi_ts, s.window = "periodic", robust = TRUE)
autoplot(stl_fit) +
  labs(title = "STL Decomposition – CPI Annual Inflation") +
  theme_minimal(base_size = 12)

adf_result  <- adf.test(cpi_ts, alternative = "stationary")
kpss_result <- kpss.test(cpi_ts, null = "Level")

results_tbl <- tibble(
  Test = c("ADF (H0: unit root)", "KPSS (H0: stationary)"),
  Statistic = c(round(adf_result$statistic, 4), round(kpss_result$statistic, 4)),
  `p-value`  = c(round(adf_result$p.value, 4), round(kpss_result$p.value, 4)),
  Conclusion = c(
    ifelse(adf_result$p.value  < 0.05, "Reject H0 → Stationary", "Fail to Reject H0 → Non-Stationary"),
    ifelse(kpss_result$p.value < 0.05, "Reject H0 → Non-Stationary", "Fail to Reject H0 → Stationary")
  )
)

kable(results_tbl, caption = "Stationarity Test Results")

# If non-stationary, check first difference
cpi_diff <- diff(cpi_ts)
adf_diff  <- adf.test(cpi_diff, alternative = "stationary")
cat("ADF on first-differenced series — p-value:", round(adf_diff$p.value, 4), "\n")

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
acf(cpi_ts,  lag.max = 20, main = "ACF – CPI y/y%")
pacf(cpi_ts, lag.max = 20, main = "PACF – CPI y/y%")
acf(cpi_diff,  lag.max = 20, main = "ACF – First Difference")
pacf(cpi_diff, lag.max = 20, main = "PACF – First Difference")

par(mfrow = c(1, 1))

fc <- forecast(arima_fit, h = 8, level = c(80, 95))

autoplot(fc) +
  geom_hline(yintercept = c(1, 3), linetype = "dashed", color = "darkgreen", linewidth = 0.8) +
  annotate("text", x = time(cpi_ts)[length(cpi_ts)] + 1,
           y = 2, label = "RBNZ Target\n1–3%", color = "darkgreen", size = 3.5) +
  labs(title = "ARIMA Forecast – CPI Annual Inflation (y/y%)",
       subtitle = paste("Model:", arima_fit$arma %>% paste(collapse = ",")),
       x = NULL, y = "y/y %",
       caption = "Shaded bands: 80% and 95% prediction intervals") +
  theme_minimal(base_size = 13)

fc_df <- data.frame(
  Quarter   = as.character(time(fc$mean)),
  Forecast  = round(fc$mean, 2),
  Lower_80  = round(fc$lower[, 1], 2),
  Upper_80  = round(fc$upper[, 1], 2),
  Lower_95  = round(fc$lower[, 2], 2),
  Upper_95  = round(fc$upper[, 2], 2)
)
kable(fc_df, caption = "8-Quarter CPI Inflation Forecast (y/y %)")

series_list <- list(
  "CPI"          = df$CPI_yy,
  "Tradable"     = df$Tradable_yy,
  "Non-Tradable" = df$NonTradable_yy,
  "Food"         = df$Food_yy,
  "Housing"      = df$Housing_yy
)

multi_results <- map_dfr(names(series_list), function(nm) {
  vals <- series_list[[nm]]
  vals <- vals[!is.na(vals)]
  ts_obj <- ts(vals, frequency = 4)
  fit <- auto.arima(ts_obj, stepwise = TRUE, approximation = TRUE)
  tibble(
    Series = nm,
    Model  = paste0("ARIMA(", paste(arimaorder(fit), collapse = ","), ")"),
    AICc   = round(fit$aicc, 2),
    RMSE   = round(sqrt(mean(residuals(fit)^2)), 4)
  )
})

kable(multi_results, caption = "Auto ARIMA Models – All Key Series")

recent <- tail(df, 8)
kable(recent %>%
  select(Date, CPI_yy, Tradable_yy, NonTradable_yy, Food_yy, Housing_yy) %>%
  mutate(Date = as.character(Date)),
  digits = 2,
  caption = "Recent 8 Quarters – Key Inflation Metrics (y/y %)")