TAR
Libraries
library(dplyr)
library(psych)
library(readxl)
library(writexl)
library(kableExtra)
# Econometrics
library(tidyverse)
library(plm) # panel models
library(sandwich) #covariance matrixes
library(lmtest) # tests
library(xtable)# latex tables
library(stargazer) # latex regression tables
library(ggpubr) # correlation test
# for this HA
library(tsDyn)
library(tseries)
library(tidyverse)
Downloading the data
time_use_dat <- read_xlsx(
'C:/Users/Popov/Documents/Studies/NES_studies/R/Time_Series/HA2/Data.xlsx')
head(time_use_dat)
## # A tibble: 6 × 2
## Date Value
## <dttm> <chr>
## 1 1948-01-01 00:00:00 3.4
## 2 1948-02-01 00:00:00 3.8
## 3 1948-03-01 00:00:00 4.0
## 4 1948-04-01 00:00:00 3.9
## 5 1948-05-01 00:00:00 3.5
## 6 1948-06-01 00:00:00 3.6
Long Difference Threshold Variable
ut <- as.integer(time_use_dat$Value)
d <- 5 # Long difference parameter
threshold_var <- diff(ut, lag = d)
Fit a Linear AR model and select lag order using BIC
max_lag <- 10
n <- length(ut) # Number of observations
bic_values <- numeric(max_lag)
for (i in 1:max_lag) {
ar_fit <- ar(ut, order.max = i, method = "mle")
sigma2 <- ar_fit$var.pred # Residual variance estimate
k <- i + 1 # Number of parameters (lags + intercept)
bic_values[i] <- n * log(sigma2) + k * log(n) # Calculate BIC
}
# Find the lag with the lowest BIC
selected_lag <- which.min(bic_values)
Fit a Threshold Autoregressive (TAR) Model
# Adjust the thDelay to be less than the selected lag
thDelay <- min(d, selected_lag - 1)
# Fit
tar_model <- setar(ut, m = selected_lag, thDelay = thDelay, trim = 0.1)
tar_summary <- summary(tar_model)
tar_summary
##
## Non linear autoregressive model
##
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
## const.L phiL.1 phiL.2
## 0.1820360 0.7958157 0.1698977
##
## High regime:
## const.H phiH.1 phiH.2
## 1.6386814 1.0908133 -0.2869785
##
## Threshold:
## -Variable: Z(t) = + (0) X(t)+ (1)X(t-1)
## -Value: 7
## Proportion of points in low regime: 90.43% High regime: 9.57%
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.390051 -0.079176 -0.010603 0.057970 10.125008
##
## Fit:
## residuals variance = 0.2677, AIC = -1201, MAPE = 4.553%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const.L 0.182036 0.069674 2.6127 0.009130 **
## phiL.1 0.795816 0.033908 23.4698 < 2.2e-16 ***
## phiL.2 0.169898 0.035059 4.8460 1.479e-06 ***
## const.H 1.638681 0.528048 3.1033 0.001973 **
## phiH.1 1.090813 0.116487 9.3642 < 2.2e-16 ***
## phiH.2 -0.286979 0.107621 -2.6666 0.007798 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (0) X(t) + (1) X(t-1)
##
## Value: 7
Confidence Interval for the Threshold
coefficients_table <- tar_summary$coefficients
if (is.matrix(coefficients_table) || is.data.frame(coefficients_table)) {
estimates <- coefficients_table[, 1] # The first column
standard_errors <- coefficients_table[, 2] # The second column
# Confidence level
z_value <- 1.96
# Calculate confidence intervals
lower_bounds <- estimates - z_value * standard_errors
upper_bounds <- estimates + z_value * standard_errors
# Combine results in a data frame
ci_results <- data.frame(
Coefficient = rownames(coefficients_table),
Estimate = estimates,
Std_Error = standard_errors,
CI_Lower = lower_bounds,
CI_Upper = upper_bounds
)
print(ci_results)
} else {
print("The coefficients_table is not in the expected format. Please check the structure of tar_summary.")
}
## Coefficient Estimate Std_Error CI_Lower CI_Upper
## const.L const.L 0.1820360 0.06967372 0.04547548 0.31859645
## phiL.1 phiL.1 0.7958157 0.03390808 0.72935583 0.86227551
## phiL.2 phiL.2 0.1698977 0.03505932 0.10118143 0.23861397
## const.H const.H 1.6386814 0.52804776 0.60370777 2.67365499
## phiH.1 phiH.1 1.0908133 0.11648745 0.86249791 1.31912872
## phiH.2 phiH.2 -0.2869785 0.10762057 -0.49791482 -0.07604219
logistic smooth transition autoregress
# Load the necessary libraries
library(tsDyn)
max_lag <- 10
bic_values <- numeric(max_lag)
# Fit AR models and choose the one with the lowest BIC
for (i in 1:max_lag) {
ar_fit <- ar(ut, order.max = i, method = "mle")
sigma2 <- ar_fit$var.pred # Residual variance estimate
k <- i + 1 # Number of parameters (lags + intercept)
bic_values[i] <- length(ut) * log(sigma2) + k * log(length(ut)) # Calculate BIC
}
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
# Find the lag order with the lowest BIC
selected_lag <- which.min(bic_values)
cat("Selected lag order by BIC:", selected_lag, "\n")
## Selected lag order by BIC: 2
# Fit the LSTAR Model
# Use the lagged difference as the transition variable
diff_ut <- diff(ut) # Compute the first difference of 'ut'
lstar_model <- lstar(ut, m = selected_lag, d = 1, thVar = lag(diff_ut, -1))
## Using maximum autoregressive order for low regime: mL = 2
## Using maximum autoregressive order for high regime: mH = 2
## Using only first 920 elements of thVar
## Performing grid search for starting values...
## Starting values fixed: gamma = 1 , th = 0 ; SSE = 249.3713
## Grid search selected lower/upper bound gamma (was: 1 100 ]).
## Might try to widen bound with arg: 'starting.control=list(gammaInt=c(1,200))'
## Optimization algorithm converged
## Optimized values fixed for regime 2 : gamma = 1.365344 , th = -2.279359 ; SSE = 245.7466
summary(lstar_model)
## Warning in vcov.lstar(object): Hessian negative-semi definite
## Warning in sqrt(diag(vc)): созданы NaN
##
## Non linear autoregressive model
##
## LSTAR model
## Coefficients:
## Low regime:
## const.L phiL.1 phiL.2
## 6.662986 1.584704 -1.296270
##
## High regime:
## const.H phiH.1 phiH.2
## -6.7875029 -0.6794527 1.3891906
##
## Smoothing parameter: gamma = 1.365
##
## Threshold
## Variable: external
## Value: -2.279
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6978810 -0.0684895 -0.0043476 0.0597943 10.1647837
##
## Fit:
## residuals variance = 0.2677, AIC = -1199, MAPE = 4.544%
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|z|)
## const.L 6.66299 NaN NaN NaN
## phiL.1 1.58470 1.67990 0.9433 0.34551
## phiL.2 -1.29627 0.90939 -1.4254 0.15404
## const.H -6.78750 NaN NaN NaN
## phiH.1 -0.67945 1.64564 -0.4129 0.67969
## phiH.2 1.38919 0.84266 1.6486 0.09923 .
## gamma 1.36534 0.21429 6.3714 1.873e-10 ***
## th -2.27936 NaN NaN NaN
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Non-linearity test of full-order LSTAR model against full-order AR model
## F = 5.9374 ; p-value = 0.0027416
##
## Threshold
## Variable: external
# Step 4: Perform a Likelihood Ratio Test (GLR Test) to Test Linearity
# Fit a linear AR model and compare with LSTAR model
# Calculate residual sum of squares for the linear model
linear_ar_model <- ar(ut, order.max = selected_lag, method = "mle")
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): возможно,
## проблема сходимости: 'optim' выдал код= 1
rss_linear <- sum(resid(linear_ar_model)^2)
# Calculate residual sum of squares for the LSTAR model
rss_lstar <- sum(lstar_model$residuals^2)
# Likelihood ratio test statistic
n <- length(ut)
likelihood_ratio_stat <- n * (log(rss_linear) - log(rss_lstar))
# Degrees of freedom: number of additional parameters in the LSTAR model
df <- 1
# Calculate the p-value
p_value <- 1 - pchisq(likelihood_ratio_stat, df)
cat("Likelihood Ratio Statistic:", likelihood_ratio_stat, "\n")
## Likelihood Ratio Statistic: NA
cat("p-value:", p_value, "\n")
## p-value: NA
# Wald Test for Transition Parameters
# Check if the transition parameter is significantly different from zero
# Extract parameter estimates from the summary
lstar_summary <- summary(lstar_model)
## Warning in vcov.lstar(object): Hessian negative-semi definite
## Warning in vcov.lstar(object): созданы NaN
coefs <- lstar_summary$coefficients
## Assuming the transition parameter is the parameter of interest
#transition_param <- coefs["gamma", "Estimate"]
#transition_se <- coefs["gamma", "Std. Error"]
## Wald test statistic for the transition parameter
#wald_stat <- (transition_param / transition_se)^2
#wald_p_value <- 1 - pchisq(wald_stat, df)
## Print Wald test results
#cat("Wald Test Statistic for Transition Parameter:", wald_stat, "\n")
#cat("p-value:", wald_p_value, "\n")