The kardl package is an R tool for estimating symmetric
and asymmetric Autoregressive Distributed Lag (ARDL) and Nonlinear ARDL
(NARDL) models, designed for econometricians and researchers analyzing
cointegration and dynamic relationships in time series data. It offers
flexible model specifications, allowing users to include deterministic
variables, asymmetric effects for short- and long-run dynamics, and
trend components. The package supports customizable lag structures,
model selection criteria (AIC, BIC, AICc, HQ), and parallel processing
for computational efficiency. Key features include:
asym(), asymL(), and asymS() to
model asymmetric effects in short- and long-run dynamics, and
deterministic() for dummy variables."quick", "grid",
"grid_custom") or user-defined lags.This vignette demonstrates how to use the kardl()
function to estimate an asymmetric ARDL model, perform diagnostic tests,
and visualize results, using economic data from Turkey.
kardl in R can easily be installed from its CRAN
repository:
Alternatively, you can use the devtools package to load
directly from GitHub:
# Install required packages
install.packages(c("stats", "msm", "lmtest", "nlWaldTest", "car", "strucchange", "utils"))
# Install kardl from GitHub
install.packages("devtools")
devtools::install_github("karamelikli/kardl")Load the package:
This example estimates an asymmetric ARDL model to analyze the
dynamics of exchange rate pass-through to domestic prices in Turkey,
using a sample dataset (imf_example_data) with variables
for Consumer Price Index (CPI), Exchange Rate (ER), Producer Price Index
(PPI), and a COVID-19 dummy variable.
Assume imf_example_data contains monthly data for CPI,
ER, PPI, and a COVID dummy variable. We prepare the data by ensuring
proper formatting and adding the dummy variable. We retrieve data from
the IMF’s International Financial Statistics (IFS) dataset and prepare
it for analysis.
Note: The imf_example_data is a placeholder for
demonstration purposes. You should replace it with your actual dataset.
The data can be loaded by readxl or other data import
functions.
We define the model formula using R’s formula syntax, incorporating
asymmetric effects and deterministic variables. We use
asym() for variables with both short- and long-run
asymmetry, deterministic() for fixed effects, and
trend for a linear time trend.
We estimate the ARDL model using different mode settings
to demonstrate flexibility in lag selection. The kardl()
function supports various modes: "grid",
"grid_custom", "quick", or a user-defined lag
vector.
mode = "grid"The "grid" mode evaluates all lag combinations up to
maxlag and provides console feedback.
# Set model options
kardl_set(criterion = "BIC", differentAsymLag = TRUE, data=imf_example_data)
# Estimate model with grid mode
kardl_model <- kardl(data=imf_example_data,model= MyFormula, maxlag = 4, mode = "grid")## ==============================================================
## KARDL Results
## The lags of the model is:
##                       AIC               BIC              AICc                HQ
## lag           2,1,0,3,0         1,1,0,3,0         1,0,0,0,0         2,1,0,3,0
## value -5.55341824428271 -5.39678297453557 -4.65456961489876 -5.49030935302091
## 
## The estimated model's formula is:
## L0.d.CPI~L1.CPI+L1.ER_POS+L1.ER_NEG+L1.PPI_POS+L1.PPI_NEG+L1.d.CPI+L0.d.ER_POS+L1.d.ER_POS+L0.d.ER_NEG+L0.d.PPI_POS+L1.d.PPI_POS+L2.d.PPI_POS+L3.d.PPI_POS+L0.d.PPI_NEG+covid+trend
## 
## The coeffients estimation's results are:
##   (Intercept)        L1.CPI     L1.ER_POS     L1.ER_NEG    L1.PPI_POS 
## -0.0386633545 -0.0121524327  0.0110491088  0.0252653399  0.0517030999 
##    L1.PPI_NEG      L1.d.CPI   L0.d.ER_POS   L1.d.ER_POS   L0.d.ER_NEG 
##  0.0451043051  0.3340367351  0.1111220281  0.0937503220 -0.0026590963 
##  L0.d.PPI_POS  L1.d.PPI_POS  L2.d.PPI_POS  L3.d.PPI_POS  L0.d.PPI_NEG 
##  0.0474101941  0.0021468066 -0.0519928228 -0.0517409420  0.0057550123 
##         covid         trend 
##  0.0033274526 -0.0002952073 
## ===============================================================
## ==============================================================
## KARDL Summary
## 
## 
## Call:
## lm(formula = as.formula(fmodel), data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.050478 -0.008129 -0.000904  0.006918  0.102836 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.0386634  0.0227251  -1.701 0.089572 .  
## L1.CPI       -0.0121524  0.0047287  -2.570 0.010494 *  
## L1.ER_POS     0.0110491  0.0051559   2.143 0.032652 *  
## L1.ER_NEG     0.0252653  0.0079538   3.177 0.001594 ** 
## L1.PPI_POS    0.0517031  0.0096244   5.372 1.25e-07 ***
## L1.PPI_NEG    0.0451043  0.0107631   4.191 3.35e-05 ***
## L1.d.CPI      0.3340367  0.0399191   8.368 7.54e-16 ***
## L0.d.ER_POS   0.1111220  0.0180412   6.159 1.63e-09 ***
## L1.d.ER_POS   0.0937503  0.0181990   5.151 3.88e-07 ***
## L0.d.ER_NEG  -0.0026591  0.0474028  -0.056 0.955291    
## L0.d.PPI_POS  0.0474102  0.0160401   2.956 0.003284 ** 
## L1.d.PPI_POS  0.0021468  0.0144158   0.149 0.881684    
## L2.d.PPI_POS -0.0519928  0.0143424  -3.625 0.000322 ***
## L3.d.PPI_POS -0.0517409  0.0137160  -3.772 0.000183 ***
## L0.d.PPI_NEG  0.0057550  0.0135155   0.426 0.670452    
## covid         0.0033275  0.0050899   0.654 0.513621    
## trend        -0.0002952  0.0002472  -1.194 0.233112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01483 on 448 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.6516, Adjusted R-squared:  0.6392 
## F-statistic: 52.38 on 16 and 448 DF,  p-value: < 2.2e-16
## 
## ===============================================================
Specify custom lags to bypass automatic lag selection:
kardl_model2 <- kardl(data=imf_example_data, MyFormula, mode = c(2, 1, 1, 3, 0))
# View results
kardl_model2$properLag##     CPI  ER_POS  ER_NEG PPI_POS PPI_NEG 
##       2       1       1       3       0
Use the . operator to include all variables except the
dependent variable:
## ==============================================================
## KARDL Results
## The lags of the model is:
##                       AIC               BIC              AICc                HQ
## lag               1,2,3             1,0,0             1,0,0             1,2,0
## value -3.30318556115875 -3.22112627740211 -2.97831477517044 -3.26494116658086
## 
## The estimated model's formula is:
## L0.d.ER~L1.ER+L1.CPI+L1.PPI+L1.d.ER+L0.d.CPI+L0.d.PPI+covid
## 
## The coeffients estimation's results are:
##  (Intercept)        L1.ER       L1.CPI       L1.PPI      L1.d.ER     L0.d.CPI 
## -0.069790957 -0.025177677  0.022530297 -0.003051024  0.078941014  0.742452528 
##     L0.d.PPI        covid 
## -0.019897958  0.020270715 
## ===============================================================
The LagCriteria component contains lag combinations and
their criterion values. We visualize these to compare model selection
criteria (AIC, BIC, HQ).
library(dplyr)
library(tidyr)
library(ggplot2)
# Convert LagCriteria to a data frame
LagCriteria <- as.data.frame(kardl_model[["LagCriteria"]])
colnames(LagCriteria) <- c("lag", "AIC", "BIC", "AICc", "HQ")
LagCriteria <- LagCriteria %>% mutate(across(c(AIC, BIC, HQ), as.numeric))
# Pivot to long format
LagCriteria_long <- LagCriteria %>%
  select(-AICc) %>%
  pivot_longer(cols = c(AIC, BIC, HQ), names_to = "Criteria", values_to = "Value")
# Find minimum values
min_values <- LagCriteria_long %>%
  group_by(Criteria) %>%
  slice_min(order_by = Value) %>%
  ungroup()
# Plot
ggplot(LagCriteria_long, aes(x = lag, y = Value, color = Criteria, group = Criteria)) +
  geom_line() +
  geom_point(data = min_values, aes(x = lag, y = Value), color = "red", size = 3, shape = 8) +
  geom_text(data = min_values, aes(x = lag, y = Value, label = lag), vjust = 1.5, color = "black", size = 3.5) +
  labs(title = "Lag Criteria Comparison", x = "Lag Configuration", y = "Criteria Value") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))We calculate long-run coefficients using
kardl_longrun(), which standardizes coefficients by
dividing them by the negative of the dependent variable’s long-run
parameter.
## ___________________________________________________
## Adjusted KARDL Results (long-run results)
##                 Estimate Std. Error   t value     Pr(>|t|)    
## (Intercept) -3.1815321  0.6750740 -4.712864 3.265356e-06 ***
## L1.ER_POS    0.9092096  0.2083414  4.364038 1.586695e-05 ***
## L1.ER_NEG    2.0790356  0.5669276  3.667197 2.746682e-04 ***
## L1.PPI_POS   4.2545473  1.7282246  2.461802 1.419976e-02   *
## L1.PPI_NEG   3.7115453  1.3351815  2.779806 5.668108e-03  **
## ___________________________________________________
##   Signif. codes: 0.001 = ***, 0.01 = **, 0.05 = *, 0.1 = .
## ___________________________________________________
The asymmetrytest() function performs Wald tests to
assess short- and long-run asymmetry in the model.
ast <- imf_example_data %>% kardl(CPI ~ ER + PPI + asym(ER + PPI) + deterministic(covid) + trend, mode = c(1, 2, 3, 0, 1)) %>% asymmetrytest()
ast## Asymmetries in the long run
## 
##             F         P
## ER  2.1649064 0.1418984
## PPI 0.7152717 0.3981528
## ________________________________
## Asymmetries in the short run
## 
##             F         P
## ER  1.6798112 0.1956198
## PPI 0.7582018 0.3843603
## $H0
## [1] "- Coef(L1.ER_POS)/Coef(L1.CPI) = - Coef(L1.ER_NEG)/Coef(L1.CPI)"  
## [2] "- Coef(L1.PPI_POS)/Coef(L1.CPI) = - Coef(L1.PPI_NEG)/Coef(L1.CPI)"
## 
## $H1
## [1] "- Coef(L1.ER_POS)/Coef(L1.CPI) ≠ - Coef(L1.ER_NEG)/Coef(L1.CPI)"  
## [2] "- Coef(L1.PPI_POS)/Coef(L1.CPI) ≠ - Coef(L1.PPI_NEG)/Coef(L1.CPI)"
## Asymmetries in the long run
## H0: - Coef(L1.ER_POS)/Coef(L1.CPI) = - Coef(L1.ER_NEG)/Coef(L1.CPI)
## H1: - Coef(L1.ER_POS)/Coef(L1.CPI) ≠ - Coef(L1.ER_NEG)/Coef(L1.CPI)
## 
## H0: - Coef(L1.PPI_POS)/Coef(L1.CPI) = - Coef(L1.PPI_NEG)/Coef(L1.CPI)
## H1: - Coef(L1.PPI_POS)/Coef(L1.CPI) ≠ - Coef(L1.PPI_NEG)/Coef(L1.CPI)
## 
##             F         P df1 df2
## ER  2.1649064 0.1418984   1 446
## PPI 0.7152717 0.3981528   1 446
## _____________________________
## 
## Asymmetries in the long run
## H0: Coef(L0.d.ER_POS) + Coef(L1.d.ER_POS) + Coef(L2.d.ER_POS) = Coef(L0.d.ER_NEG) + Coef(L1.d.ER_NEG) + Coef(L2.d.ER_NEG) + Coef(L3.d.ER_NEG)
## H1: Coef(L0.d.ER_POS) + Coef(L1.d.ER_POS) + Coef(L2.d.ER_POS) ≠ Coef(L0.d.ER_NEG) + Coef(L1.d.ER_NEG) + Coef(L2.d.ER_NEG) + Coef(L3.d.ER_NEG)
## 
## H0: Coef(L0.d.PPI_POS) = Coef(L0.d.PPI_NEG) + Coef(L1.d.PPI_NEG)
## H1: Coef(L0.d.PPI_POS) ≠ Coef(L0.d.PPI_NEG) + Coef(L1.d.PPI_NEG)
## 
##             F         P DF    Sum.of.Sq ResDF1 ResDF2      RSS1      RSS2
## ER  1.6798112 0.1956198  1 0.0003860722    447    446 0.1028906 0.1025045
## PPI 0.7582018 0.3843603  1 0.0001742581    447    446 0.1026788 0.1025045
We perform cointegration tests to assess long-term relationships
using pssf(), psst(), narayan(),
banerjee(), and recmt().
The pssf() function tests for cointegration using the
Pesaran, Shin, and Smith F Bound test.
A <- kardl_model %>% pssf(case = 3, signif_level = "0.05")
cat(paste0("The F statistic = ", A$statistic, " where k = ", A$k, "."))## The F statistic = 11.2705611215356 where k = 2.
## 
## We found 'Cointegration' at 0.05.
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   2.00   4.19   5.06   4.87   5.85   5.49   6.59   6.34   7.52
## Method: PesaranF. Case: 5
## H0: Coef(L1.CPI) = Coef(L1.ER_POS) = Coef(L1.ER_NEG) = Coef(L1.PPI_POS) = Coef(L1.PPI_NEG) = 0
## H1: Coef(L1.CPI) ≠ Coef(L1.ER_POS) ≠ Coef(L1.ER_NEG) ≠ Coef(L1.PPI_POS) ≠ Coef(L1.PPI_NEG)≠ 0
## The F statistics = 11.2705611215356 where k = 2. 
## The proboblity is = 0.0000000003
## The results: There is Cointegration
## Significance level: 0.05
## ___________________________________________________
## Critical values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   2.00   4.19   5.06   4.87   5.85   5.49   6.59   6.34   7.52
The psst() function tests the significance of the lagged
dependent variable’s coefficient.
A <- kardl_model %>% psst(case = 3, signif_level = "0.05")
cat(paste0("The t statistic = ", A$statistic, " where k = ", A$k, "."))## The t statistic = -2.56993172541683 where k = 4.
## 
## We found 'No Cointegration' at 0.05.
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   3.00  -3.13  -3.84  -3.41  -4.16  -3.65  -4.42  -3.96  -4.73
## Method: Pesarant. Case: 5
## H0: Coef(L1.CPI) = 0
## H1: Coef(L1.CPI)≠ 0
## The t statistics = -2.56993172541683 where k = 4.
## The results: There is No Cointegration
## Significance level: 0.05
## ___________________________________________________
## Critical values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   3.00  -3.13  -3.84  -3.41  -4.16  -3.65  -4.42  -3.96  -4.73
The narayan() function is tailored for small sample
sizes.
A <- kardl_model %>% narayan(case = 3, signif_level = "0.05")
cat(paste0("The F statistic = ", A$statistic, " where k = ", A$k, "."))## The F statistic = 11.2705611215356 where k = 2.
## 
## We found 'Cointegration' at 0.05.
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##  2.000  4.307  5.223  5.067  6.103     NA     NA  6.730  8.053
## Method: Narayan. Case: 5
## H0: Coef(L1.CPI) = Coef(L1.ER_POS) = Coef(L1.ER_NEG) = Coef(L1.PPI_POS) = Coef(L1.PPI_NEG) = 0
## H1: Coef(L1.CPI) ≠ Coef(L1.ER_POS) ≠ Coef(L1.ER_NEG) ≠ Coef(L1.PPI_POS) ≠ Coef(L1.PPI_NEG)≠ 0
## The F statistics = 11.2705611215356 where k = 2.
## The proboblity is = 0.0000000003
## The results: There is Cointegration
## Significance level: 0.05
## ___________________________________________________
## Critical values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##  2.000  4.307  5.223  5.067  6.103     NA     NA  6.730  8.053
The banerjee() function is designed for small datasets
(≤100 observations).
A <- kardl_model %>% banerjee(signif_level = "0.05")
cat(paste0("The ECM parameter = ", A$coef, ", k = ", A$k, ", t statistic = ", A$statistic, "."))## The ECM parameter = -0.00867941225158909, k = 4, t statistic = -2.07931701669587.
## 
## We found 'No Cointegration' at 0.05.
##  0.01  0.05  0.10  0.25 
## -5.07 -4.38 -4.02 -3.46
## Method: Banerjee. Case: NULL
## H0: Coef(L1.CPI) = 0
## H1: Coef(L1.CPI)≠ 0
## The coeffient of ECM is -0.00867941225158909 , t=-2.07931701669587 k=4
## The results: There is No Cointegration
## Significance level: 0.05
## ___________________________________________________
## Critical values are as follows:
##  0.01  0.05  0.10  0.25 
## -5.07 -4.38 -4.02 -3.46
The recmt() function tests for cointegration using an
Error Correction Model.
## 
## RESM test. Case: 5
## The t statistics = 3.73402572442152 where k = 4.
## 
## Warnings:
## Warning: Trend is used in the model. The case was set to 5.
## Method: recmt. Case: 5
## H0: Coef(EcmRes) = 0
## H1: Coef(EcmRes)≠ 0
## Attention: This function was performed in two steps.
## In the first step, the ECM was calculated and in the second step, the PSS test was performed.
## 
## The coeffient of ECM is 0.00261085240978681 , t=3.73402572442152 k=4
## The results: There is No Cointegration
## Significance level: 
## ___________________________________________________
## Critical values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   3.00  -3.13  -3.84  -3.41  -4.16  -3.65  -4.42  -3.96  -4.73
The archtest() function checks for autoregressive
conditional heteroskedasticity in the model residuals.
## ARCH test
## 
## F = 8.372629, p-value = 0.0002681842, df1 = 2, df2 = 460
## 
## 
## Call:
## lm(formula = as.formula(paste0("resid~", paste("resid", 1:q, 
##     sep = "", collapse = "+"))), data = ndata)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0006394 -0.0001738 -0.0001317 -0.0000137  0.0103827 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.747e-04  3.091e-05   5.652 2.79e-08 ***
## resid1       1.902e-01  4.660e-02   4.081 5.28e-05 ***
## resid2      -2.159e-02  4.649e-02  -0.464    0.643    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0006081 on 460 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.03512,    Adjusted R-squared:  0.03093 
## F-statistic: 8.373 on 2 and 460 DF,  p-value: 0.0002682
We demonstrate how to customize prefixes and suffixes for asymmetric
variables using kardl_set().
# Set custom prefixes and suffixes
kardl_reset()
kardl_set(AsymPrefix = c("asyP_", "asyN_"), AsymSuffix = c("_PP", "_NN"))
kardl_custom <- kardl(data=imf_example_data, MyFormula, mode = "grid_custom")
kardl_custom$properLag##         CPI  asyP_ER_PP  asyN_ER_NN asyP_PPI_PP asyN_PPI_NN 
##           2           1           0           3           0
kardl(data, model, maxlag, mode, ...):
data: A time series dataset (e.g., a data frame with
CPI, ER, PPI).model: A formula specifying the long-run equation,
e.g.,
y ~ x + z + asym(z) + asymL(x2 + x3) + asymS(x3 + x4) + deterministic(dummy1 + dummy2) + trend.
Supports:
asym(): Asymmetric effects for both short- and long-run
dynamics.asymL(): Long-run asymmetric variables.asymS(): Short-run asymmetric variables.deterministic(): Fixed dummy variables.trend: Linear time trend.maxlag: Maximum number of lags (default: 4). Use
smaller values (e.g., 2) for small datasets, larger values (e.g., 8) for
long-term dependencies.mode: Estimation mode:
"quick": Verbose output for interactive use."grid": Verbose output with lag optimization."grid_custom": Silent, efficient execution.c(1, 2, 4, 5) or
c(CPI = 2, ER_POS = 3, ER_NEG = 1, PPI = 3)).inputs,
finalModel, start_time, end_time,
properLag, TimeSpan, OptLag,
LagCriteria, type (“kardlmodel”).kardl_set(...): Configures options
like criterion (AIC, BIC, AICc, HQ),
differentAsymLag, AsymPrefix,
AsymSuffix, ShortCoef, and
LongCoef. Use kardl_get() to retrieve settings
and kardl_reset() to restore defaults.
kardl_longrun(model): Calculates
standardized long-run coefficients, returning type
(“kardl_longrun”), coef, delta_se,
results, and starsDesc.
asymmetrytest(model): Performs Wald
tests for short- and long-run asymmetry, returning
Lhypotheses, Lwald, Shypotheses,
Swald, and type (“asymmetrytest”).
pssf(model, case, signif_level):
Performs the Pesaran, Shin, and Smith F Bound test for cointegration,
supporting cases 1–5 and significance levels (“auto”, 0.01, 0.025, 0.05,
0.1, 0.10).
psst(model, case, signif_level):
Performs the PSS t Bound test, focusing on the lagged dependent
variable’s coefficient.
narayan(model, case, signif_level):
Conducts the Narayan test for cointegration, optimized for small samples
(cases 2–5).
banerjee(model, signif_level):
Performs the Banerjee cointegration test for small datasets (≤100
observations).
recmt(data, model, maxlag, mode, case, signif_level, ...):
Conducts the Restricted ECM test for cointegration, with similar
parameters to kardl() and case/significance level
options.
archtest(resid, q): Tests for ARCH
effects in model residuals, returning type,
statistic, parameter, p.value,
and Fval.
For detailed documentation, use ?kardl,
?kardl_set, ?kardl_longrun,
?asymmetrytest, ?pssf, ?psst,
?narayan, ?banerjee, ?recmt, or
?archtest.
The kardl package is a versatile tool for econometric
analysis, offering robust support for symmetric and asymmetric
ARDL/NARDL modeling, cointegration tests, stability diagnostics, and
heteroskedasticity checks. Its flexible formula specification, lag
optimization, and support for parallel processing make it ideal for
studying complex economic relationships. For more information, visit https://github.com/karamelikli/kardl
or contact the authors at hakperest@gmail.com.