library(kardl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
library(kardl)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
# Read the pi data (assumed to be the same for all sections)
# Initialize an empty data frame to store results
results <- data.frame(
Section = character(),
Intercept = numeric(),
L1.lpiger_POS = numeric(),
L1.lpiger_NEG = numeric(),
L1.lrer_POS = numeric(),
L1.lrer_NEG = numeric(),
F_stat = numeric(),
RESET_F = numeric(),
PSS_F = numeric(),
LM_stat = numeric(),
Wald_lpiger = numeric(),
Wald_lrer = numeric(),
R_squared = numeric(),
# Add columns for p-values
p_Intercept = numeric(),
p_L1.lpiger_POS = numeric(),
p_L1.lpiger_NEG = numeric(),
p_L1.lrer_POS = numeric(),
p_L1.lrer_NEG = numeric(),
p_F_stat = numeric(),
p_RESET_F = numeric(),
p_PSS_F = numeric(),
p_LM_stat = numeric(),
p_Wald_lpiger = numeric(),
p_Wald_lrer = numeric(),
stringsAsFactors = FALSE
)
for (i in 1:21) {
tryCatch({
cat("Processing section", i, "\n")
section_data <- read_excel(paste0("section", i, ".xlsx"))
df <- bind_cols(section_data, pi)
# Filter out invalid rows
df <- df %>%
filter(EXPTR > 0 & !is.na(EXPTR) &
RER > 0 & !is.na(RER) &
TurPI > 0 & !is.na(TurPI) &
GerPI > 0 & !is.na(GerPI)) %>%
mutate(exptr = ifelse(EXPTR < 1000, EXPTR * 1000, EXPTR))
# Transformations
df$lexptr <- log(df$exptr)
df$lrer <- log(df$RER)
df$lpiger <- log(df$TurPI)
df$lpitr <- log(df$GerPI)
# Check logged variables
if (any(is.na(df$lexptr)) || any(is.na(df$lrer)) || any(is.na(df$lpiger)) ||
any(is.infinite(df$lexptr)) || any(is.infinite(df$lrer)) || any(is.infinite(df$lpiger))) {
warning(paste("Section", i, "still has NA/Inf in model variables after filtering"))
next
}
# Fit the model
kardl_m <- kardl(df, lexptr ~ trend + asym(lpiger + lrer), maxlag = 3, mode = "performance")
# Extract coefficients
coef_names <- c("(Intercept)", "L1.lpiger_POS", "L1.lpiger_NEG", "L1.lrer_POS", "L1.lrer_NEG")
coefs <- kardl_m$finalModel$model$coefficients[coef_names]
# Extract p-values for coefficients
sum_model <- summary(kardl_m$finalModel$model)
# Make sure to handle potential missing coefficients
coef_p_values <- rep(NA, length(coef_names))
names(coef_p_values) <- coef_names
for (name in coef_names) {
if (name %in% rownames(sum_model$coefficients)) {
coef_p_values[name] <- sum_model$coefficients[name, "Pr(>|t|)"]
}
}
# Get F-statistic, p-value, and R-squared
F_stat <- sum_model$fstatistic[1]
F_p_value <- pf(sum_model$fstatistic[1], sum_model$fstatistic[2], sum_model$fstatistic[3], lower.tail = FALSE)
R_squared <- sum_model$r.squared
# Perform RESET test
reset_test <- resettest(kardl_m$finalModel$model)
RESET_F <- reset_test$statistic[1]
RESET_p <- reset_test$p.value
# Perform PSS F-test for cointegration
hypothesis <- c("L1.lexptr = 0", "L1.lpiger_POS = 0", "L1.lpiger_NEG = 0",
"L1.lrer_POS = 0", "L1.lrer_NEG = 0")
wald_test <- linearHypothesis(kardl_m$finalModel$model, hypothesis)
PSS_F <- wald_test[2, "F"]
PSS_p <- wald_test[2, "Pr(>F)"]
# Perform LM test for serial correlation
lm_test <- bgtest(kardl_m$finalModel$model, order = 1)
LM_stat <- lm_test$statistic
LM_p <- lm_test$p.value
# Perform asymmetry tests and extract Wald F-statistics and p-values
ast <- asymmetrytest(kardl_m)
# The asymmetry test results have columns "F" and "P" (not "Pr(>F)")
Wald_lpiger <- ast$Lwald["lpiger", "F"]
Wald_lpiger_p <- ast$Lwald["lpiger", "P"]
Wald_lrer <- ast$Lwald["lrer", "F"]
Wald_lrer_p <- ast$Lwald["lrer", "P"]
# Collect results into a data frame row
row <- data.frame(
Section = paste("Section", i),
Intercept = coefs["(Intercept)"],
L1.lpiger_POS = coefs["L1.lpiger_POS"],
L1.lpiger_NEG = coefs["L1.lpiger_NEG"],
L1.lrer_POS = coefs["L1.lrer_POS"],
L1.lrer_NEG = coefs["L1.lrer_NEG"],
F_stat = F_stat,
RESET_F = RESET_F,
PSS_F = PSS_F,
LM_stat = LM_stat,
Wald_lpiger = Wald_lpiger,
Wald_lrer = Wald_lrer,
R_squared = R_squared,
# Add p-values
p_Intercept = coef_p_values["(Intercept)"],
p_L1.lpiger_POS = coef_p_values["L1.lpiger_POS"],
p_L1.lpiger_NEG = coef_p_values["L1.lpiger_NEG"],
p_L1.lrer_POS = coef_p_values["L1.lrer_POS"],
p_L1.lrer_NEG = coef_p_values["L1.lrer_NEG"],
p_F_stat = F_p_value,
p_RESET_F = RESET_p,
p_PSS_F = PSS_p,
p_LM_stat = LM_p,
p_Wald_lpiger = Wald_lpiger_p,
p_Wald_lrer = Wald_lrer_p,
stringsAsFactors = FALSE
)
# Append to results
results <- rbind(results, row)
cat("Completed section", i, "\n")
}, error = function(e) {
warning(paste("Error processing section", i, ":", e$message))
})
}
## Processing section 1
## Warning in value[[3L]](cond): Error processing section 1 : `path` does not
## exist: 'section1.xlsx'
## Processing section 2
## Warning in value[[3L]](cond): Error processing section 2 : `path` does not
## exist: 'section2.xlsx'
## Processing section 3
## Warning in value[[3L]](cond): Error processing section 3 : `path` does not
## exist: 'section3.xlsx'
## Processing section 4
## Warning in value[[3L]](cond): Error processing section 4 : `path` does not
## exist: 'section4.xlsx'
## Processing section 5
## Warning in value[[3L]](cond): Error processing section 5 : `path` does not
## exist: 'section5.xlsx'
## Processing section 6
## Warning in value[[3L]](cond): Error processing section 6 : `path` does not
## exist: 'section6.xlsx'
## Processing section 7
## Warning in value[[3L]](cond): Error processing section 7 : `path` does not
## exist: 'section7.xlsx'
## Processing section 8
## Warning in value[[3L]](cond): Error processing section 8 : `path` does not
## exist: 'section8.xlsx'
## Processing section 9
## Warning in value[[3L]](cond): Error processing section 9 : `path` does not
## exist: 'section9.xlsx'
## Processing section 10
## Warning in value[[3L]](cond): Error processing section 10 : `path` does not
## exist: 'section10.xlsx'
## Processing section 11
## Warning in value[[3L]](cond): Error processing section 11 : `path` does not
## exist: 'section11.xlsx'
## Processing section 12
## Warning in value[[3L]](cond): Error processing section 12 : `path` does not
## exist: 'section12.xlsx'
## Processing section 13
## Warning in value[[3L]](cond): Error processing section 13 : `path` does not
## exist: 'section13.xlsx'
## Processing section 14
## Warning in value[[3L]](cond): Error processing section 14 : `path` does not
## exist: 'section14.xlsx'
## Processing section 15
## Warning in value[[3L]](cond): Error processing section 15 : `path` does not
## exist: 'section15.xlsx'
## Processing section 16
## Warning in value[[3L]](cond): Error processing section 16 : `path` does not
## exist: 'section16.xlsx'
## Processing section 17
## Warning in value[[3L]](cond): Error processing section 17 : `path` does not
## exist: 'section17.xlsx'
## Processing section 18
## Warning in value[[3L]](cond): Error processing section 18 : `path` does not
## exist: 'section18.xlsx'
## Processing section 19
## Warning in value[[3L]](cond): Error processing section 19 : `path` does not
## exist: 'section19.xlsx'
## Processing section 20
## Warning in value[[3L]](cond): Error processing section 20 : `path` does not
## exist: 'section20.xlsx'
## Processing section 21
## Warning in value[[3L]](cond): Error processing section 21 : `path` does not
## exist: 'section21.xlsx'
# Convert numeric columns to numeric (in case of any coercion)
numeric_cols <- names(results)[2:23]
results[, numeric_cols] <- lapply(results[, numeric_cols], as.numeric)
# Create a function to add significance stars
add_stars <- function(value, p_value) {
if (is.na(value) || is.na(p_value)) return(as.character(round(value, 4)))
stars <- ""
if (p_value < 0.01) stars <- "***"
else if (p_value < 0.05) stars <- "**"
else if (p_value < 0.1) stars <- "*"
return(paste0(round(value, 4), stars))
}
# Create a new results table with stars
results_with_stars <- data.frame(
Section = results$Section,
Intercept = mapply(add_stars, results$Intercept, results$p_Intercept),
L1.lpiger_POS = mapply(add_stars, results$L1.lpiger_POS, results$p_L1.lpiger_POS),
L1.lpiger_NEG = mapply(add_stars, results$L1.lpiger_NEG, results$p_L1.lpiger_NEG),
L1.lrer_POS = mapply(add_stars, results$L1.lrer_POS, results$p_L1.lrer_POS),
L1.lrer_NEG = mapply(add_stars, results$L1.lrer_NEG, results$p_L1.lrer_NEG),
F_stat = mapply(add_stars, results$F_stat, results$p_F_stat),
RESET_F = mapply(add_stars, results$RESET_F, results$p_RESET_F),
PSS_F = mapply(add_stars, results$PSS_F, results$p_PSS_F),
LM_stat = mapply(add_stars, results$LM_stat, results$p_LM_stat),
Wald_lpiger = mapply(add_stars, results$Wald_lpiger, results$p_Wald_lpiger),
Wald_lrer = mapply(add_stars, results$Wald_lrer, results$p_Wald_lrer),
R_squared = round(results$R_squared, 4)
)
# View the resulting table with significance indicators
print(results_with_stars)
## [1] Section R_squared
## <0 rows> (or 0-length row.names)
# Optional: Write to CSV
write.csv(results_with_stars, "kardl_results_with_significance.csv", row.names = FALSE)