Nothing
R/ardlBound.R
In dLagM: Time Series Regression Models with Distributed Lag Models
Defines functions
ardlBound
appendList
cusumSq
Documented in
ardlBound
#' @import stats
#' @import graphics
#' @importFrom lmtest dwtest
#' @importFrom lmtest bptest
#' @importFrom lmtest waldtest
#' @importFrom strucchange efp
#' @importFrom strucchange recresid
#' @importFrom grDevices dev.new
#' @importFrom grDevices graphics.off
#' @importFrom utils tail
#' @importFrom sandwich NeweyWest
cusumSq <- function(model, alpha = 0.05){
# The following table is Table 1 of
# Tests for Serial Correlation in Regression Analysis Based on the Periodogram of LeastSquares Residuals
# Author(s): J. Durbin
# Source: Biometrika, Vol. 56, No. 1 (Mar., 1969), pp. 1-15
# Published by: Oxford University Press on behalf of Biometrika Trust
# Stable URL: https://www.jstor.org/stable/2334686
# Accessed: 13-08-2019 01:44 UTC
table <- matrix( c(NA, 0.1 , 0.05, 0.025, 0.01, 0.005,
1 , 0.4 , 0.45 , 0.475 , 0.49 , 0.495 ,
2 , 0.35044 , 0.44306 , 0.50855 , 0.56667 , 0.59596 ,
3 , 0.35477 , 0.41811 , 0.46702 , 53456 , 0.579 ,
4 , 0.33435 , 0.39075 , 0.44641 , 0.50495 , 0.5421 ,
5 , 0.31556 , 0.37359 , 0.42174 , 0.47692 , 0.51576 ,
6 , 0.30244 , 0.35522 , 0.40045 , 0.4544 , 0.48988 ,
7 , 0.28991 , 33905 , 0.38294 , 0.43337 , 0.46761 ,
8 , 0.27828 , 0.32538 , 0.36697 , 0.41522 , 0.44819 ,
9 , 0.26794 , 0.31325 , 0.35277 , 0.39922 , 0.43071 ,
10 , 0.25884 , 0.30221 , 0.34022 , 0.38481 , 0.41517 ,
11 , 0.25071 , 0.29227 , 0.32894 , 0.37187 , 0.40122 ,
12 , 0.24325 , 0.2833 , 0.31869 , 0.36019 , 0.38856 ,
13 , 0.23639 , 0.27515 , 0.30935 , 0.34954 , 0.37703 ,
14 , 0.2301 , 0.26767 , 0.30081 , 0.3398 , 0.36649 ,
15 , 0.2243 , 0.26077 , 0.29296 , 0.33083 , 0.35679 ,
16 , 0.21895 , 0.25439 , 0.2857 , 0.32256 , 0.34784 ,
17 , 0.21397 , 0.24847 , 0.27897 , 0.31489 , 0.33953 ,
18 , 0.20933 , 0.24296 , 0.2727 , 0.30775 , 0.33181 ,
19 , 0.20498 , 23781 , 0.26685 , 0.30108 , 0.32459 ,
20 , 0.20089 , 0.23298 , 0.26137 , 0.29484 , 0.31784 ,
21 , 0.19705 , 0.22844 , 0.25622 , 0.28898 , 0.31149 ,
22 , 0.19343 , 0.22416 , 0.25136 , 0.28346 , 0.30552 ,
23 , 0.19001 , 0.22012 , 0.24679 , 0.27825 , 0.29989 ,
24 , 18677 , 0.2163 , 0.24245 , 0.27333 , 0.29456 ,
25 , 18370 , 0.21268 , 0.23835 , 0.26866 , 0.28951 ,
26 , 0.18077 , 0.20924 , 0.23445 , 0.26423 , 0.28472 ,
27 , 0.17799 , 20596 , 0.23074 , 0.26001 , 0.28016 ,
28 , 0.17533 , 0.20283 , 0.22721 , 0.256 , 0.27582 ,
29 , 0.1728 , 19985 , 0.22383 , 0.25217 , 0.27168 ,
30 , 0.17037 , 0.197 , 0.22061 , 0.24851 , 0.26772 ,
31 , 0.16805 , 0.19427 , 0.21752 , 0.24501 , 0.26393 ,
32 , 0.16582 , 19166 , 0.21457 , 0.24165 , 0.2603 ,
33 , 0.16368 , 0.18915 , 0.21173 , 0.23843 , 0.25683 ,
34 , 0.16162 , 0.18674 , 0.20901 , 0.23534 , 0.25348 ,
35 , 0.15964 , 0.18442 , 0.20639 , 0.23237 , 0.25027 ,
36 , 0.15774 , 18218 , 0.20387 , 0.22951 , 0.24718 ,
37 , 0.1559 , 0.18003 , 0.20144 , 0.22676 , 0.24421 ,
38 , 0.15413 , 0.17796 , 0.1991 , 0.2241 , 0.24134 ,
39 , 0.15242 , 17595 , 0.19684 , 0.22154 , 0.23857 ,
40 , 0.15076 , 0.17402 , 0.19465 , 0.21906 , 0.23589 ,
42 , 0.14761 , 0.17034 , 0.1905 , 0.21436 , 0.23081 ,
43 , 0.14611 , 0.16858 , 0.18852 , 0.21212 , 0.22839 ,
44 , 0.14466 , 0.16688 , 0.18661 , 0.20995 , 0.22605 ,
45 , 0.14325 , 0.16524 , 0.18475 , 20785 , 0.22377 ,
46 , 0.14188 , 0.16364 , 0.18295 , 0.20581 , 0.22157 ,
47 , 0.14055 , 0.16208 , 0.1812 , 0.20383 , 0.21943 ,
48 , 0.13926 , 0.16058 , 0.1795 , 0.2019 , 0.21735 ,
49 , 0.138 , 0.15911 , 0.17785 , 0.20003 , 0.21534 ,
50 , 13678 , 0.15769 , 0.17624 , 0.19822 , 0.21337 ,
51 , 0.13559 , 0.1563 , 0.17468 , 0.19645 , 0.21146 ,
52 , 0.13443 , 0.15495 , 0.17316 , 0.19473 , 0.20961 ,
53 , 0.1333 , 0.15363 , 0.17168 , 0.19305 , 0.2078 ,
54 , 0.13221 , 0.15235 , 0.17024 , 0.19142 , 0.20604 ,
55 , 0.13113 , 0.1511 , 0.16884 , 0.18983 , 0.20432 ,
56 , 0.13009 , 0.14989 , 0.16746 , 0.18828 , 0.20265 ,
57 , 0.12907 , 0.1487 , 0.16613 , 0.18677 , 0.20101 ,
58 , 0.12807 , 0.14754 , 0.16482 , 0.18529 , 0.19942 ,
59 , 0.1271 , 0.14641 , 0.16355 , 0.18385 , 0.19786 ,
60 , 0.12615 , 0.1453 , 0.1623 , 0.18245 , 0.19635 ,
62 , 0.12431 , 0.14316 , 0.1599 , 0.17973 , 0.19341 ,
64 , 0.12255 , 0.14112 , 0.1576 , 0.17713 , 0.19061 ,
66 , 0.12087 , 0.13916 , 0.1554 , 0.17464 , 0.18792 ,
68 , 0.11926 , 0.13728 , 0.15329 , 0.17226 , 0.18535 ,
70 , 0.11771 , 0.13548 , 0.15127 , 0.16997 , 0.18288 ,
72 , 0.11622 , 0.13375 , 0.14932 , 0.16777 , 0.18051 ,
74 , 0.11479 , 0.13208 , 0.14745 , 0.16566 , 0.17823 ,
76 , 0.11341 , 0.13048 , 0.14565 , 0.16363 , 0.17604 ,
78 , 0.11208 , 0.12894 , 0.14392 , 16167 , 0.17392 ,
80 , 0.11079 , 0.12745 , 0.14224 , 0.15978 , 0.17188 ,
82 , 0.10955 , 0.12601 , 0.14063 , 0.15795 , 0.16992 ,
84 , 0.10835 , 0.12462 , 0.13907 , 0.15619 , 0.16802 ,
86 , 0.10719 , 0.12327 , 0.13756 , 0.15449 , 0.16618 ,
88 , 10607 , 0.12197 , 0.1361 , 15284 , 0.1644 ,
90 , 0.10499 , 0.12071 , 0.13468 , 0.15124 , 0.16268 ,
92 , 0.10393 , 0.11949 , 0.13331 , 0.1497 , 0.16101 ,
94 , 0.10291 , 0.11831 , 0.13198 , 0.1482 , 0.594 ,
96 , 0.10192 , 0.11716 , 0.1307 , 0.14674 , 0.15783 ,
98 , 0.10096 , 0.11604 , 0.12944 , 0.14533 , 0.15631 ,
100 , 0.10002 , 0.11496 , 0.12823 , 0.14396 , 0.15483 ), nrow = 80, ncol = 6, byrow = TRUE)
# model = model1.1$ECM$EC.model
Wr <- recresid(model, data = model$model, type = "Rec-CUSUM")
k <- start(Wr)[1] - 1 #length(model$coefficients)
t <- end(Wr)[1] + 1
if (((t-k) %% 2) != 0){
n1 <- min( 0.5 * (t-k) - 1.5, 100)
n2 <- min( 0.5 * (t-k) - 0.5, 100)
c01 <- table[which((table[,1] == n1)), which(table[1,] == alpha/2)]
c02 <- table[which((table[,1] == n2)), which(table[1,] == alpha/2)]
c0 <- mean(c(c01,c02))
} else {
n <- min ( 0.5 * (t-k)-1, 100)
c0 <- table[which((table[,1] == n)), which(table[1,] == alpha/2)]
}
WrSq <- array(NA, length(Wr))
lb <- array(NA, length(Wr))
ub <- array(NA, length(Wr))
for (i in 1:length(Wr)){
WrSq[i] <- sum(Wr[1:i]^2)/sum(Wr^2)
lb[i] <- -c0 + (start(Wr)[1] + i - 1 - k)/(t-k)
ub[i] <- c0 + (start(Wr)[1] + i - 1 - k)/(t-k)
}
# Remove
cex = 1.1
par(family = "Arial")
par(cex.lab=cex, cex.main = cex, cex.axis = cex)
# Remove
plot(ts(WrSq, start = start(Wr)), ylim = c(min(lb)-0.1,max(ub)+0.1), ylab="CUSUM of squared residuals", main ="Recursive CUSUM of squares test")
lines(ts(lb, start = start(Wr)), col = "red")
lines(ts(ub, start = start(Wr)), col = "red")
abline(h=0, col = rgb(0,0,0), lty = 1, lwd=1.05)
}
appendList <- function(list1, list2){
for (i in names(list1)[which(names(list1) == names(list2))]){
list1[[i]] = sort(c(list1[[i]], list2[[i]]))
}
for (i in names(list2)[which(names(list2) != names(list1))]){
list1[[i]] = list2[[i]]
}
return(list1)
}
#' @export
ardlBound <- function(data = NULL , formula = NULL , case = 3 , p = NULL , remove = NULL, autoOrder = FALSE ,
HAC = FALSE,
ic = c("AIC", "BIC", "MASE", "GMRAE") ,
max.p = 15, max.q = 15, ECM = TRUE, stability = TRUE){
if (is.null(data)) stop("Enter data by data argument.")
if (is.null(formula)) stop("A formula object showing the dependent and indepdenent series must be entered.")
vars <- all.vars(formula)
NumVar <- length(vars)
if ( NumVar < ncol(data)){
data <- data[,vars]
}
if (is.null(p) | (autoOrder == TRUE) ){
cat(" ", "\n")
cat("Orders being calculated with max.p =", max.p , "and max.q =", max.q, "...\n\n")
orders <- ardlBoundOrders(data = data , formula = formula, ic = "AIC" , max.p = max.p, max.q = max.q )
cat("Autoregressive order:" , orders$q + 1, "and p-orders:" , unlist(orders$p) + 1 , "\n")
cat("------------------------------------------------------", "\n")
p <- data.frame(orders$q , orders$p ) + 1
}
if (!is.data.frame(p)){
p <- array(p , NumVar)
p <- data.frame(t(p))
colnames(p) <- vars
} else {
colnames(p) <- vars
if (length(p) < NumVar) stop("Enter an integer or and array of lenght equal to the number of variables in the formula!")
}
caseType <- switch(case, "no intercept, no trend", "restricted intercert, no trend (not supported)",
"unrestricted intercert, no trend", "unrestricted intercept, restricted trend (not supported)",
"unrestricted intercept, unrestricted trend")
diffData <- apply(data , 2 , diff)
colnames(diffData) <- paste0("d" , colnames(diffData))
data <- cbind(data[2:nrow(data),],diffData)
if (max(p[2:NumVar]) == 1){
max.p <- 2
} else {
max.p <- max(p[2:NumVar])
}
rem.p <- list()
for (i in 1:NumVar){
if (max.p > 2){
rem.p[[vars[i] ]] <- c(0,2:(max.p-1))
}else {
rem.p[[vars[i] ]] <- c(0)
}
}
if (sum((p[2:NumVar] - max.p) == 0) == (NumVar-1)){
removeP <- NULL
if (ECM == TRUE) {
remP <- list()
if ((max.p-1) >= 2){
remP[["ec"]] <- c(0, 2:(max.p-1))
} else if ((max.p-1) == 1){
remP[["ec"]] <- c(0)
}
removeP = list(p = remP )
}
} else {
remP <- list()
reduce <- which((p[2:NumVar] - max.p) != 0) + 1
for ( j in reduce){
remP[[paste0("d" ,vars[j])]] <- c((max.p-1):p[[j]])
rem.p[[paste0("d" ,vars[j])]] <- c((max.p-1):p[[j]])
}
if (ECM == TRUE) {
if ((max.p-1) >= 2){
remP[["ec"]] <- c(0, 2:(max.p-1))
} else if ((max.p-1) == 1){
remP[["ec"]] <- c(0)
}
}
# if (((max.p-1) >= 2) & (ECM == TRUE)) remP[["ec"]] <- c(0, 2:(max.p-1))
removeP = list(p = remP )
}
removeP$p <- appendList(removeP$p, remove$p)
removeP$q <- remove$q
rem.p <- appendList(rem.p, remove$p)
removeP2 <- list(p = rem.p, q = remove$q )
formula1 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(vars ,
paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (case == 1){
formula1 <- update(formula1, ~. -1) # Case 1 requires this
formula2 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
# print(formula1)
# print(formula2)
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) ,
p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[2:NumVar] / modelFull$model$coefficients[1])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 2){
formula1 <- update(formula1, ~. +1)
formula2 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) ,
p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar])) + (modelFull$model$coefficients[1] / modelFull$model$coefficients[2])
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 3){
formula2 <- as.formula(paste0("d" , vars[1] , " ~ +1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ +1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 4){
trend <- c(1:nrow(data))
data <- cbind(data , trend)
formula1 <- as.formula(paste0(formula1 , " + trend"))
formula2 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
rem.p[["trend"]] <- c(1:(max.p-1))
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar])) + modelFull$model$coefficients["trend.t"] * trend / modelFull$model$coefficients[2]
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 5){
trend <- c(1:nrow(data))
data <- cbind(data , trend)
formula1 <- as.formula(paste0(formula1 , " + trend"))
formula2 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") , " + trend" ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
removeP$p[["trend"]] <- c(1:(max.p-1))
removeP2$p[["trend"]] <- c(1:(max.p-1))
rem.p[["trend"]] <- c(1:(max.p-1))
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2 )
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") , " + trend" ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
}
bg <- NULL
bg <- lmtest::bgtest(modelFull$model, type = "F")
cat("\n Breusch-Godfrey Test for the autocorrelation in residuals:\n")
print(bg)
if (!is.na(bg$p.value)){
if (bg$p.value < 0.05){
cat("The p-value of Breusch-Godfrey test for the autocorrelation in residuals: ", bg$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
res <- modelFull$model$residuals
cat("\n Ljung-Box Test for the autocorrelation in residuals:\n")
lb <- Box.test(res, type = c("Ljung-Box"))
print(lb)
if (!is.na(lb$p.value)){
if (lb$p.value < 0.05){
cat("The p-value of Ljung-Box test for the autocorrelation in residuals: ", lb$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
bp <- NULL
bp <- lmtest::bptest(modelFull$model)
cat("\n Breusch-Pagan Test for the homoskedasticity of residuals:\n")
print(bp)
if (!is.na(bp$p.value)){
if (bp$p.value < 0.05){
cat("The p-value of Breusch-Pagan test for the homoskedasticity of residuals: ", bp$p.value , "< 0.05!\n" )
# stop("Bounds test procedure is stopped due to significant heteroskedasticity in residuals.")
}
}
cat("------------------------------------------------------", "\n")
sp <- NULL
sp <- shapiro.test(modelFull$model$residual)
cat("\n Shapiro-Wilk test of normality of residuals:\n")
print(sp)
if (!is.na(sp$p.value)){
if (sp$p.value < 0.05){
cat("The p-value of Shapiro-Wilk test normality of residuals: ", sp$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
tryCatch(
{if (HAC == TRUE) {
cat("Newey-West HAC covariance matrix estimators are used for testing.", "\n")
neweyCM <- NeweyWest(modelFull$model)
Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model, vcov = neweyCM)$F[2]
} else {
Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model)$F[2]
}
},
error = function(e) {
Fvalue <<- anova( modelNull$model , modelFull$model)$F[2]
}
)
# Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model)$F[2]
# if (is.null(k)) k = (NumVar - 1)
k = (NumVar - 1)
pssbounds(obs = nrow(data) , fstat = Fvalue , case = case, k = k )
if ((stability == TRUE) & (ECM == TRUE)){
cat("------------------------------------------------------", "\n")
reset <- NULL
reset <- lmtest::resettest(modelECM$model)
cat("\n Ramsey's RESET Test for model specification:\n")
print(reset)
if (reset$p.value < 0.05){
cat("the p-value of RESET test: ", reset$p.value , "< 0.05!\n" )
}
cat("------------------------------------------------------", "\n")
cusum.test <- efp(modelECM$model, data = modelECM$model$model, type = "Rec-CUSUM")
mosum.test <- efp(modelECM$model, data = modelECM$model$model, type = "Rec-MOSUM")
graphics.off()
plot.new()
# Remove
cex = 1.1
par(family = "Arial")
par(cex.lab=cex, cex.main = cex, cex.axis = cex)
# Remove
par(oma = c(4, 1, 1, 1))
two <- FALSE
if (sum(c(any(is.nan(cusum.test$process)), any(is.nan(mosum.test$process)))) >0){
par(mfrow = c(1,2))
two <- TRUE
} else {
par(mfrow = c(2,2))
}
if (!any(is.nan(cusum.test$process))){
plot(cusum.test)
}
cusumSq(modelECM$model)
if(!any(is.nan(mosum.test$process))){
plot(mosum.test)
}
if (two){
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot.new()
legend("bottom", c("Recursive residuals", "5% limits"), xpd = TRUE, horiz = TRUE, inset = 0, bty = "n",
col = c("black" , "red"), lwd=2, cex = 1)
} else {
# plot(0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
plot.new()
legend("top", c("Recursive residuals", "5% limits"), col = c("black" , "red"),
xpd = TRUE, bty = "o", lwd = 2, cex = 1)
}
}
if (ECM){
cat("------------------------------------------------------", "\n")
cat("Error Correction Model Output: \n")
summary(modelECM)
cat("------------------------------------------------------", "\n")
cat("Long-run coefficients: \n")
if (case == 1){
print(modelFull$model$coefficients[1:NumVar])
} else {
print(modelFull$model$coefficients[2:(NumVar+1)])
}
cat(" ", "\n")
# text <- paste0("EC = " , names(modelFull$model$coefficients[2]), " - (")
# for ( i in 3:(NumVar+1)){
# text <- paste0(text, "(", modelFull$model$coefficients[i]/modelFull$model$coefficients[2], " * ",
# names(modelFull$model$coefficients[i]) , ") + ")
# }
# cat(text, "\n")
}
if (ECM){
return(list(model = list(modelNull = modelNull, modelFull = modelFull), F.stat = Fvalue , p = p , k = k, bg = bg, lb = lb , bp = bp, sp = sp,
ECM = list(EC.t = ec, EC.model = modelECM$model , EC.beta = ecm.beta, EC.data = data.ecm), ARDL.model = modelFull$model ) )
} else {
return(list(model = list(modelNull = modelNull, modelFull = modelFull), ARDL.model = modelFull$model, F.stat = Fvalue , p = p , k = k, bg = bg, lb = lb , bp = bp, sp = sp ) )
}
}
Try the dLagM package in your browser
Any scripts or data that you put into this service are public.
dLagM documentation built on Oct. 2, 2023, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ardlBound appendList cusumSq
Documented in ardlBound
#' @import stats
#' @import graphics
#' @importFrom lmtest dwtest
#' @importFrom lmtest bptest
#' @importFrom lmtest waldtest
#' @importFrom strucchange efp
#' @importFrom strucchange recresid
#' @importFrom grDevices dev.new
#' @importFrom grDevices graphics.off
#' @importFrom utils tail
#' @importFrom sandwich NeweyWest
cusumSq <- function(model, alpha = 0.05){
# The following table is Table 1 of
# Tests for Serial Correlation in Regression Analysis Based on the Periodogram of LeastSquares Residuals
# Author(s): J. Durbin
# Source: Biometrika, Vol. 56, No. 1 (Mar., 1969), pp. 1-15
# Published by: Oxford University Press on behalf of Biometrika Trust
# Stable URL: https://www.jstor.org/stable/2334686
# Accessed: 13-08-2019 01:44 UTC
table <- matrix( c(NA, 0.1 , 0.05, 0.025, 0.01, 0.005,
1 , 0.4 , 0.45 , 0.475 , 0.49 , 0.495 ,
2 , 0.35044 , 0.44306 , 0.50855 , 0.56667 , 0.59596 ,
3 , 0.35477 , 0.41811 , 0.46702 , 53456 , 0.579 ,
4 , 0.33435 , 0.39075 , 0.44641 , 0.50495 , 0.5421 ,
5 , 0.31556 , 0.37359 , 0.42174 , 0.47692 , 0.51576 ,
6 , 0.30244 , 0.35522 , 0.40045 , 0.4544 , 0.48988 ,
7 , 0.28991 , 33905 , 0.38294 , 0.43337 , 0.46761 ,
8 , 0.27828 , 0.32538 , 0.36697 , 0.41522 , 0.44819 ,
9 , 0.26794 , 0.31325 , 0.35277 , 0.39922 , 0.43071 ,
10 , 0.25884 , 0.30221 , 0.34022 , 0.38481 , 0.41517 ,
11 , 0.25071 , 0.29227 , 0.32894 , 0.37187 , 0.40122 ,
12 , 0.24325 , 0.2833 , 0.31869 , 0.36019 , 0.38856 ,
13 , 0.23639 , 0.27515 , 0.30935 , 0.34954 , 0.37703 ,
14 , 0.2301 , 0.26767 , 0.30081 , 0.3398 , 0.36649 ,
15 , 0.2243 , 0.26077 , 0.29296 , 0.33083 , 0.35679 ,
16 , 0.21895 , 0.25439 , 0.2857 , 0.32256 , 0.34784 ,
17 , 0.21397 , 0.24847 , 0.27897 , 0.31489 , 0.33953 ,
18 , 0.20933 , 0.24296 , 0.2727 , 0.30775 , 0.33181 ,
19 , 0.20498 , 23781 , 0.26685 , 0.30108 , 0.32459 ,
20 , 0.20089 , 0.23298 , 0.26137 , 0.29484 , 0.31784 ,
21 , 0.19705 , 0.22844 , 0.25622 , 0.28898 , 0.31149 ,
22 , 0.19343 , 0.22416 , 0.25136 , 0.28346 , 0.30552 ,
23 , 0.19001 , 0.22012 , 0.24679 , 0.27825 , 0.29989 ,
24 , 18677 , 0.2163 , 0.24245 , 0.27333 , 0.29456 ,
25 , 18370 , 0.21268 , 0.23835 , 0.26866 , 0.28951 ,
26 , 0.18077 , 0.20924 , 0.23445 , 0.26423 , 0.28472 ,
27 , 0.17799 , 20596 , 0.23074 , 0.26001 , 0.28016 ,
28 , 0.17533 , 0.20283 , 0.22721 , 0.256 , 0.27582 ,
29 , 0.1728 , 19985 , 0.22383 , 0.25217 , 0.27168 ,
30 , 0.17037 , 0.197 , 0.22061 , 0.24851 , 0.26772 ,
31 , 0.16805 , 0.19427 , 0.21752 , 0.24501 , 0.26393 ,
32 , 0.16582 , 19166 , 0.21457 , 0.24165 , 0.2603 ,
33 , 0.16368 , 0.18915 , 0.21173 , 0.23843 , 0.25683 ,
34 , 0.16162 , 0.18674 , 0.20901 , 0.23534 , 0.25348 ,
35 , 0.15964 , 0.18442 , 0.20639 , 0.23237 , 0.25027 ,
36 , 0.15774 , 18218 , 0.20387 , 0.22951 , 0.24718 ,
37 , 0.1559 , 0.18003 , 0.20144 , 0.22676 , 0.24421 ,
38 , 0.15413 , 0.17796 , 0.1991 , 0.2241 , 0.24134 ,
39 , 0.15242 , 17595 , 0.19684 , 0.22154 , 0.23857 ,
40 , 0.15076 , 0.17402 , 0.19465 , 0.21906 , 0.23589 ,
42 , 0.14761 , 0.17034 , 0.1905 , 0.21436 , 0.23081 ,
43 , 0.14611 , 0.16858 , 0.18852 , 0.21212 , 0.22839 ,
44 , 0.14466 , 0.16688 , 0.18661 , 0.20995 , 0.22605 ,
45 , 0.14325 , 0.16524 , 0.18475 , 20785 , 0.22377 ,
46 , 0.14188 , 0.16364 , 0.18295 , 0.20581 , 0.22157 ,
47 , 0.14055 , 0.16208 , 0.1812 , 0.20383 , 0.21943 ,
48 , 0.13926 , 0.16058 , 0.1795 , 0.2019 , 0.21735 ,
49 , 0.138 , 0.15911 , 0.17785 , 0.20003 , 0.21534 ,
50 , 13678 , 0.15769 , 0.17624 , 0.19822 , 0.21337 ,
51 , 0.13559 , 0.1563 , 0.17468 , 0.19645 , 0.21146 ,
52 , 0.13443 , 0.15495 , 0.17316 , 0.19473 , 0.20961 ,
53 , 0.1333 , 0.15363 , 0.17168 , 0.19305 , 0.2078 ,
54 , 0.13221 , 0.15235 , 0.17024 , 0.19142 , 0.20604 ,
55 , 0.13113 , 0.1511 , 0.16884 , 0.18983 , 0.20432 ,
56 , 0.13009 , 0.14989 , 0.16746 , 0.18828 , 0.20265 ,
57 , 0.12907 , 0.1487 , 0.16613 , 0.18677 , 0.20101 ,
58 , 0.12807 , 0.14754 , 0.16482 , 0.18529 , 0.19942 ,
59 , 0.1271 , 0.14641 , 0.16355 , 0.18385 , 0.19786 ,
60 , 0.12615 , 0.1453 , 0.1623 , 0.18245 , 0.19635 ,
62 , 0.12431 , 0.14316 , 0.1599 , 0.17973 , 0.19341 ,
64 , 0.12255 , 0.14112 , 0.1576 , 0.17713 , 0.19061 ,
66 , 0.12087 , 0.13916 , 0.1554 , 0.17464 , 0.18792 ,
68 , 0.11926 , 0.13728 , 0.15329 , 0.17226 , 0.18535 ,
70 , 0.11771 , 0.13548 , 0.15127 , 0.16997 , 0.18288 ,
72 , 0.11622 , 0.13375 , 0.14932 , 0.16777 , 0.18051 ,
74 , 0.11479 , 0.13208 , 0.14745 , 0.16566 , 0.17823 ,
76 , 0.11341 , 0.13048 , 0.14565 , 0.16363 , 0.17604 ,
78 , 0.11208 , 0.12894 , 0.14392 , 16167 , 0.17392 ,
80 , 0.11079 , 0.12745 , 0.14224 , 0.15978 , 0.17188 ,
82 , 0.10955 , 0.12601 , 0.14063 , 0.15795 , 0.16992 ,
84 , 0.10835 , 0.12462 , 0.13907 , 0.15619 , 0.16802 ,
86 , 0.10719 , 0.12327 , 0.13756 , 0.15449 , 0.16618 ,
88 , 10607 , 0.12197 , 0.1361 , 15284 , 0.1644 ,
90 , 0.10499 , 0.12071 , 0.13468 , 0.15124 , 0.16268 ,
92 , 0.10393 , 0.11949 , 0.13331 , 0.1497 , 0.16101 ,
94 , 0.10291 , 0.11831 , 0.13198 , 0.1482 , 0.594 ,
96 , 0.10192 , 0.11716 , 0.1307 , 0.14674 , 0.15783 ,
98 , 0.10096 , 0.11604 , 0.12944 , 0.14533 , 0.15631 ,
100 , 0.10002 , 0.11496 , 0.12823 , 0.14396 , 0.15483 ), nrow = 80, ncol = 6, byrow = TRUE)
# model = model1.1$ECM$EC.model
Wr <- recresid(model, data = model$model, type = "Rec-CUSUM")
k <- start(Wr)[1] - 1 #length(model$coefficients)
t <- end(Wr)[1] + 1
if (((t-k) %% 2) != 0){
n1 <- min( 0.5 * (t-k) - 1.5, 100)
n2 <- min( 0.5 * (t-k) - 0.5, 100)
c01 <- table[which((table[,1] == n1)), which(table[1,] == alpha/2)]
c02 <- table[which((table[,1] == n2)), which(table[1,] == alpha/2)]
c0 <- mean(c(c01,c02))
} else {
n <- min ( 0.5 * (t-k)-1, 100)
c0 <- table[which((table[,1] == n)), which(table[1,] == alpha/2)]
}
WrSq <- array(NA, length(Wr))
lb <- array(NA, length(Wr))
ub <- array(NA, length(Wr))
for (i in 1:length(Wr)){
WrSq[i] <- sum(Wr[1:i]^2)/sum(Wr^2)
lb[i] <- -c0 + (start(Wr)[1] + i - 1 - k)/(t-k)
ub[i] <- c0 + (start(Wr)[1] + i - 1 - k)/(t-k)
}
# Remove
cex = 1.1
par(family = "Arial")
par(cex.lab=cex, cex.main = cex, cex.axis = cex)
# Remove
plot(ts(WrSq, start = start(Wr)), ylim = c(min(lb)-0.1,max(ub)+0.1), ylab="CUSUM of squared residuals", main ="Recursive CUSUM of squares test")
lines(ts(lb, start = start(Wr)), col = "red")
lines(ts(ub, start = start(Wr)), col = "red")
abline(h=0, col = rgb(0,0,0), lty = 1, lwd=1.05)
}
appendList <- function(list1, list2){
for (i in names(list1)[which(names(list1) == names(list2))]){
list1[[i]] = sort(c(list1[[i]], list2[[i]]))
}
for (i in names(list2)[which(names(list2) != names(list1))]){
list1[[i]] = list2[[i]]
}
return(list1)
}
#' @export
ardlBound <- function(data = NULL , formula = NULL , case = 3 , p = NULL , remove = NULL, autoOrder = FALSE ,
HAC = FALSE,
ic = c("AIC", "BIC", "MASE", "GMRAE") ,
max.p = 15, max.q = 15, ECM = TRUE, stability = TRUE){
if (is.null(data)) stop("Enter data by data argument.")
if (is.null(formula)) stop("A formula object showing the dependent and indepdenent series must be entered.")
vars <- all.vars(formula)
NumVar <- length(vars)
if ( NumVar < ncol(data)){
data <- data[,vars]
}
if (is.null(p) | (autoOrder == TRUE) ){
cat(" ", "\n")
cat("Orders being calculated with max.p =", max.p , "and max.q =", max.q, "...\n\n")
orders <- ardlBoundOrders(data = data , formula = formula, ic = "AIC" , max.p = max.p, max.q = max.q )
cat("Autoregressive order:" , orders$q + 1, "and p-orders:" , unlist(orders$p) + 1 , "\n")
cat("------------------------------------------------------", "\n")
p <- data.frame(orders$q , orders$p ) + 1
}
if (!is.data.frame(p)){
p <- array(p , NumVar)
p <- data.frame(t(p))
colnames(p) <- vars
} else {
colnames(p) <- vars
if (length(p) < NumVar) stop("Enter an integer or and array of lenght equal to the number of variables in the formula!")
}
caseType <- switch(case, "no intercept, no trend", "restricted intercert, no trend (not supported)",
"unrestricted intercert, no trend", "unrestricted intercept, restricted trend (not supported)",
"unrestricted intercept, unrestricted trend")
diffData <- apply(data , 2 , diff)
colnames(diffData) <- paste0("d" , colnames(diffData))
data <- cbind(data[2:nrow(data),],diffData)
if (max(p[2:NumVar]) == 1){
max.p <- 2
} else {
max.p <- max(p[2:NumVar])
}
rem.p <- list()
for (i in 1:NumVar){
if (max.p > 2){
rem.p[[vars[i] ]] <- c(0,2:(max.p-1))
}else {
rem.p[[vars[i] ]] <- c(0)
}
}
if (sum((p[2:NumVar] - max.p) == 0) == (NumVar-1)){
removeP <- NULL
if (ECM == TRUE) {
remP <- list()
if ((max.p-1) >= 2){
remP[["ec"]] <- c(0, 2:(max.p-1))
} else if ((max.p-1) == 1){
remP[["ec"]] <- c(0)
}
removeP = list(p = remP )
}
} else {
remP <- list()
reduce <- which((p[2:NumVar] - max.p) != 0) + 1
for ( j in reduce){
remP[[paste0("d" ,vars[j])]] <- c((max.p-1):p[[j]])
rem.p[[paste0("d" ,vars[j])]] <- c((max.p-1):p[[j]])
}
if (ECM == TRUE) {
if ((max.p-1) >= 2){
remP[["ec"]] <- c(0, 2:(max.p-1))
} else if ((max.p-1) == 1){
remP[["ec"]] <- c(0)
}
}
# if (((max.p-1) >= 2) & (ECM == TRUE)) remP[["ec"]] <- c(0, 2:(max.p-1))
removeP = list(p = remP )
}
removeP$p <- appendList(removeP$p, remove$p)
removeP$q <- remove$q
rem.p <- appendList(rem.p, remove$p)
removeP2 <- list(p = rem.p, q = remove$q )
formula1 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(vars ,
paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (case == 1){
formula1 <- update(formula1, ~. -1) # Case 1 requires this
formula2 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
# print(formula1)
# print(formula2)
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) ,
p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[2:NumVar] / modelFull$model$coefficients[1])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 2){
formula1 <- update(formula1, ~. +1)
formula2 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) ,
p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar])) + (modelFull$model$coefficients[1] / modelFull$model$coefficients[2])
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ - 1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 3){
formula2 <- as.formula(paste0("d" , vars[1] , " ~ +1 + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ +1 + ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 4){
trend <- c(1:nrow(data))
data <- cbind(data , trend)
formula1 <- as.formula(paste0(formula1 , " + trend"))
formula2 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
rem.p[["trend"]] <- c(1:(max.p-1))
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2)
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar])) + modelFull$model$coefficients["trend.t"] * trend / modelFull$model$coefficients[2]
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
} else if (case == 5){
trend <- c(1:nrow(data))
data <- cbind(data , trend)
formula1 <- as.formula(paste0(formula1 , " + trend"))
formula2 <- as.formula(paste0("d" , vars[1] , " ~ " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") , " + trend" ))
if (p[[vars[1]]] == 1){
p[[vars[1]]] <- 2
removeP <- appendList(removeP, list( q = c(1)))
removeP2$q <- c(1)
}
removeP$p[["trend"]] <- c(1:(max.p-1))
removeP2$p[["trend"]] <- c(1:(max.p-1))
rem.p[["trend"]] <- c(1:(max.p-1))
modelFull <- ardlDlm(formula = formula1, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP2 )
modelNull <- ardlDlm(formula = formula2, data = data.frame(data) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
if (length(modelFull$model$residuals) != length(modelNull$model$residuals)){
strt <- nrow(data) - length(modelNull$model$residuals)
modelNull <- ardlDlm(formula = formula2, data = data.frame(tail(data,(length(modelFull$model$residuals) + strt))) , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
}
if ( ECM == TRUE){
mult <- as.vector(modelFull$model$coefficients[3:(NumVar+1)] / modelFull$model$coefficients[2])
ec <- data[,1] + as.vector(mult %*% t(data[,2:NumVar]))
data.ecm <- data.frame(cbind(data, ec))
formula3 <- as.formula(paste0("d" , vars[1] , " ~ ec + " , paste(c(paste0("d" , vars[2:NumVar] )) , collapse=" + ") , " + trend" ))
modelECM <- ardlDlm(formula = formula3, data = data.ecm , p = (max.p-1) , q = (p[[vars[1]]]-1) , remove = removeP )
ecm.beta <- modelECM$model$coefficients["ec.1"]
}
}
bg <- NULL
bg <- lmtest::bgtest(modelFull$model, type = "F")
cat("\n Breusch-Godfrey Test for the autocorrelation in residuals:\n")
print(bg)
if (!is.na(bg$p.value)){
if (bg$p.value < 0.05){
cat("The p-value of Breusch-Godfrey test for the autocorrelation in residuals: ", bg$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
res <- modelFull$model$residuals
cat("\n Ljung-Box Test for the autocorrelation in residuals:\n")
lb <- Box.test(res, type = c("Ljung-Box"))
print(lb)
if (!is.na(lb$p.value)){
if (lb$p.value < 0.05){
cat("The p-value of Ljung-Box test for the autocorrelation in residuals: ", lb$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
bp <- NULL
bp <- lmtest::bptest(modelFull$model)
cat("\n Breusch-Pagan Test for the homoskedasticity of residuals:\n")
print(bp)
if (!is.na(bp$p.value)){
if (bp$p.value < 0.05){
cat("The p-value of Breusch-Pagan test for the homoskedasticity of residuals: ", bp$p.value , "< 0.05!\n" )
# stop("Bounds test procedure is stopped due to significant heteroskedasticity in residuals.")
}
}
cat("------------------------------------------------------", "\n")
sp <- NULL
sp <- shapiro.test(modelFull$model$residual)
cat("\n Shapiro-Wilk test of normality of residuals:\n")
print(sp)
if (!is.na(sp$p.value)){
if (sp$p.value < 0.05){
cat("The p-value of Shapiro-Wilk test normality of residuals: ", sp$p.value , "< 0.05!\n" )
}
}
cat("------------------------------------------------------", "\n")
tryCatch(
{if (HAC == TRUE) {
cat("Newey-West HAC covariance matrix estimators are used for testing.", "\n")
neweyCM <- NeweyWest(modelFull$model)
Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model, vcov = neweyCM)$F[2]
} else {
Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model)$F[2]
}
},
error = function(e) {
Fvalue <<- anova( modelNull$model , modelFull$model)$F[2]
}
)
# Fvalue <- lmtest::waldtest( modelNull$model , modelFull$model)$F[2]
# if (is.null(k)) k = (NumVar - 1)
k = (NumVar - 1)
pssbounds(obs = nrow(data) , fstat = Fvalue , case = case, k = k )
if ((stability == TRUE) & (ECM == TRUE)){
cat("------------------------------------------------------", "\n")
reset <- NULL
reset <- lmtest::resettest(modelECM$model)
cat("\n Ramsey's RESET Test for model specification:\n")
print(reset)
if (reset$p.value < 0.05){
cat("the p-value of RESET test: ", reset$p.value , "< 0.05!\n" )
}
cat("------------------------------------------------------", "\n")
cusum.test <- efp(modelECM$model, data = modelECM$model$model, type = "Rec-CUSUM")
mosum.test <- efp(modelECM$model, data = modelECM$model$model, type = "Rec-MOSUM")
graphics.off()
plot.new()
# Remove
cex = 1.1
par(family = "Arial")
par(cex.lab=cex, cex.main = cex, cex.axis = cex)
# Remove
par(oma = c(4, 1, 1, 1))
two <- FALSE
if (sum(c(any(is.nan(cusum.test$process)), any(is.nan(mosum.test$process)))) >0){
par(mfrow = c(1,2))
two <- TRUE
} else {
par(mfrow = c(2,2))
}
if (!any(is.nan(cusum.test$process))){
plot(cusum.test)
}
cusumSq(modelECM$model)
if(!any(is.nan(mosum.test$process))){
plot(mosum.test)
}
if (two){
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot.new()
legend("bottom", c("Recursive residuals", "5% limits"), xpd = TRUE, horiz = TRUE, inset = 0, bty = "n",
col = c("black" , "red"), lwd=2, cex = 1)
} else {
# plot(0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
plot.new()
legend("top", c("Recursive residuals", "5% limits"), col = c("black" , "red"),
xpd = TRUE, bty = "o", lwd = 2, cex = 1)
}
}
if (ECM){
cat("------------------------------------------------------", "\n")
cat("Error Correction Model Output: \n")
summary(modelECM)
cat("------------------------------------------------------", "\n")
cat("Long-run coefficients: \n")
if (case == 1){
print(modelFull$model$coefficients[1:NumVar])
} else {
print(modelFull$model$coefficients[2:(NumVar+1)])
}
cat(" ", "\n")
# text <- paste0("EC = " , names(modelFull$model$coefficients[2]), " - (")
# for ( i in 3:(NumVar+1)){
# text <- paste0(text, "(", modelFull$model$coefficients[i]/modelFull$model$coefficients[2], " * ",
# names(modelFull$model$coefficients[i]) , ") + ")
# }
# cat(text, "\n")
}
if (ECM){
return(list(model = list(modelNull = modelNull, modelFull = modelFull), F.stat = Fvalue , p = p , k = k, bg = bg, lb = lb , bp = bp, sp = sp,
ECM = list(EC.t = ec, EC.model = modelECM$model , EC.beta = ecm.beta, EC.data = data.ecm), ARDL.model = modelFull$model ) )
} else {
return(list(model = list(modelNull = modelNull, modelFull = modelFull), ARDL.model = modelFull$model, F.stat = Fvalue , p = p , k = k, bg = bg, lb = lb , bp = bp, sp = sp ) )
}
}
Try the dLagM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.