Projeto Cimento
Thomas Daher
2023-11-17
Instalando os pacotes
Importando dados e deflacionando
Análise de raiz unitária para series em nível e em primeira diferença
Podemos notar que fora a série calcário que já é estacionária em nível as demais séries apresentam raiz unitária, ao calcular a primeira diferença todas se apresentam estacionárias.
Sazonalidade das séries
A presença de tendência e sazonalidade pode ser observada na decomposição da série em todas as séries. Consequentemente, é aconselhável empregar o comando auto.sarima para verificar a determinação do modelo na explicação da série.
SARIMA
Curiosamente o auto.arima determina primeira diferença para o calcario que já é estacionário em primeira ordem como observado pelo teste ADF. O metodo de Bai_perron será empregado para determinar a existência de quebra estrutural.
Zivot-Andrews test
#cimento
ts_cimento <- ts(ts_cimento, start = c(1979, 12), frequency = 12)
summary(ur.za((ts_cimento), model = c("both"), lag = NULL))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -140.829 -4.421 -0.528 4.643 86.979
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.25975 3.33100 -3.080 0.00224 **
## y.l1 0.87524 0.02193 39.918 < 2e-16 ***
## trend 0.39044 0.07220 5.408 1.20e-07 ***
## du -23.71999 4.38255 -5.412 1.17e-07 ***
## dt -0.34571 0.07041 -4.910 1.41e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.4 on 343 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9569, Adjusted R-squared: 0.9564
## F-statistic: 1902 on 4 and 343 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -5.6901
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 121#calcario
ts_calcario <- ts(ts_calcario, start = c(1979, 12), frequency = 12)
summary(ur.za((ts_calcario), model = c("both"), lag = NULL))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.372 -1.976 -0.362 1.927 36.203
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.03281 2.04169 5.894 9.04e-09 ***
## y.l1 0.78174 0.03285 23.798 < 2e-16 ***
## trend 0.11637 0.02261 5.147 4.47e-07 ***
## du -4.85195 1.40760 -3.447 0.000637 ***
## dt -0.07513 0.01955 -3.842 0.000145 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.736 on 343 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9404, Adjusted R-squared: 0.9397
## F-statistic: 1354 on 4 and 343 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -6.6443
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 123#diesel
ts_diesel <- ts(ts_diesel, start = c(1979, 12), frequency = 12)
summary(ur.za((ts_diesel), model = c("both"), lag = NULL))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -91.664 -5.037 -0.957 3.839 70.550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.69641 6.59503 4.200 3.41e-05 ***
## y.l1 0.86496 0.02239 38.628 < 2e-16 ***
## trend 0.33119 0.11724 2.825 0.00501 **
## du -31.98726 6.00292 -5.329 1.80e-07 ***
## dt -0.32191 0.11833 -2.720 0.00685 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.5 on 343 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9401, Adjusted R-squared: 0.9394
## F-statistic: 1345 on 4 and 343 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -6.0307
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 65#combustivel
ts_combustivel <- ts(ts_combustivel, start = c(1979, 12), frequency = 12)
summary(ur.za((ts_combustivel), model = c("both"), lag = NULL))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -63.578 -4.475 -0.975 3.483 47.451
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.51981 4.94734 4.350 1.80e-05 ***
## y.l1 0.89373 0.02338 38.234 < 2e-16 ***
## trend -0.04991 0.01506 -3.314 0.001018 **
## du 11.72512 2.95163 3.972 8.67e-05 ***
## dt 0.19311 0.05816 3.320 0.000996 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.74 on 343 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.979, Adjusted R-squared: 0.9788
## F-statistic: 4006 on 4 and 343 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -4.5464
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 231O cimento é estacionário com um nível de confiança de 1% e apresenta uma quebra estrutural em dezembro de 1989. Essa data específica está dentro do prazo indicado pelo CADE como o início do cartel, conforme evidenciado por documentos de reuniões realizadas entre 1987 e 1989 entre as cimenteiras. O calcário, por outro lado, também permanece estacionário, mas sofre uma ruptura em fevereiro de 1990, dois meses após a ocorrência no cimento. Essa discrepância pode sugerir que a quebra foi atribuída aos custos de produção e não ao próprio cartel. O diesel, da mesma forma, permanece estacionário e sofre uma quebra em abril de 1985, aparentemente desconectado da fratura no cimento. Ao contrário, a série envolvendo combustível não é estacionária com um nível de confiança de 10% e se rompe em fevereiro de 1999, aparentemente sem relação com a quebra no cimento. Em conclusão, embora a quebra no cimento esteja dentro do período em análise, ela poderia ser explicada pela quebra do calcário, um dos componentes utilizados em sua produção. Uma investigação mais meticulosa poderia fornecer informações sobre se uma pode realmente ser explicada pela outra, e isso deve ser prosseguido. Pois apesar de ambas as quebras ocorrerem, a proximidade das quebras em termos de tempo não implica necessariamente uma relação causal.
Analise da quebra
# Create a data frame for ggplot
df <- data.frame(date = as.Date(time(ts_cimento)), ts_cimento = as.vector(ts_cimento), ts_calcario = as.vector(ts_calcario))
# Subset the data to the desired time period (1987 to 1990)
df_subset <- subset(df, date >= as.Date("1987-01-01") & date <= as.Date("1990-12-31"))
ggplot(df_subset, aes(x = date)) +
geom_line(aes(y = ts_cimento, color = "ts_cimento"), size = 1) +
geom_line(aes(y = ts_calcario, color = "ts_calcario"), size = 1) +
labs(title = "Time Series Plot", x = "Date", y = "Value") +
scale_color_manual(name = "Series", values = c("ts_cimento" = "blue", "ts_calcario" = "red")) +
theme_minimal() +
geom_vline(xintercept = as.Date("1989-12-01"), linetype = "dashed", color = "blue") +
geom_vline(xintercept = as.Date("1990-02-01"), linetype = "dashed", color = "red")
Eu executo esse esquema com a intenção de atingir dois objetivos: Verificar se as quebras ocorrem na mesma orientação e determinar o grau de proporcionalidade. Deduzo que, apesar da ocorrência de intervalos de tempo comparáveis, não há evidências aparentes que sugiram que uma quebra cause a outra. Consequentemente, deve existir um fator exógeno adicional que contribua para esse fenômeno. Outra questão fundamental é que a quebra é negativa, o que não corresponde as ações esperadas do cartel divulgadas pelo CADE.
Lee-Strazicich function
ur.ls <- function(y, model = c("crash", "break"), breaks = 1, lags = NULL, method = c("GTOS","Fixed"), pn = 0.1, print.results = c("print", "silent")){
#Starttime
starttime <- Sys.time()
#Check sanity of the function call
if (any(is.na(y)))
stop("\nNAs in y.\n")
y <- as.vector(y)
if(pn >= 1 || pn <= 0){
stop("\n pn has to be between 0 and 1.")
}
if(method == "Fixed" && is.null(lags) == TRUE){
stop("\n If fixed lag length should be estimated, the number
\n of lags to be included should be defined explicitely.")
}
#Add lagmatrix function
lagmatrix <- function(x,max.lag){
embed(c(rep(NA,max.lag),x),max.lag+1)
}
#Add diffmatrix function
diffmatrix <- function(x,max.diff,max.lag){
#Add if condition to make it possible to differentiate between matrix and vector
if(is.vector(x) == TRUE ){
embed(c(rep(NA,max.lag),diff(x,max.lag,max.diff)),max.diff)
}
else if(is.matrix(x) == TRUE){
rbind(matrix(rep(NA,max.lag), ncol = ncol(x)), matrix(diff(x,max.lag,max.diff), ncol = ncol(x)))
}
#if matrix the result should be 0, if only a vector it should be 1
else if(as.integer(is.null(ncol(x))) == 0 ){
rbind(matrix(rep(NA,max.lag), ncol = ncol(x)), matrix(diff(x,max.lag,max.diff), ncol = ncol(x)))
}
}
#Number of observations
n <- length(y)
model <- match.arg(model)
lags <- as.integer(lags)
method <- match.arg(method)
breaks <- as.integer(breaks)
#Percentage to eliminate endpoints in the lag calculation
pn <- pn
#Critical Values for the one break test
model.one.crash.cval <- matrix(c(-4.239, -3.566, -3.211)
, nrow = 1, ncol = 3, byrow = TRUE)
colnames(model.one.crash.cval) <- c("1%","5%","10%")
model.one.break.cval <- matrix(c(.1 , -5.11, -4.5, -4.21,
.2 , -5.07, -4.47, -4.20,
.3 , -5.15, -4.45, -4.18,
.4 , -5.05, -4.50, -4.18,
.5 , -5.11, -4.51, -4.17), nrow = 5, ncol = 4, byrow = TRUE)
colnames(model.one.break.cval) <- c("lambda","1%","5%","10%")
# All critical values were derived in samples of size T = 100. Critical values
# in Model C (intercept and trend break) depend (somewhat) on the location of the
# break (λ = T_B /T) and are symmetric around λ and (1-λ). Model C critical values
# at additional break locations can be interpolated.
#
# Critical values for the endogenous two break test
# Model A - "crash" model
# invariant to the location of the crash
model.two.crash.cval <- matrix(c("LM_tau", -4.545, -3.842, -3.504,
"LM_rho", -35.726, -26.894, -22.892), nrow = 2, ncol = 4, byrow = TRUE )
colnames(model.two.crash.cval) <- c("Statistic","1%","5%","10%")
# Model C (i) - "break" model, breaks in the data generating process
# Model C (i) - "break" model invariant to the location of the crash
model.two.break.dgp.cval <- matrix(c("LM_tau", -5.823, -5.286, -4.989,
"LM_rho", -52.550, -45.531, -41.663), nrow = 2, ncol = 4, byrow = TRUE )
# Model C (ii) - "break" model, breaks are not considered in the data generating process
# Model C (ii) - "break" model depends on the location of the crash
## highest level of list is the location of the second breakpoint - so the share inside
## the matrix refers to the first breakpoint
model.two.break.tau.cval <- matrix(c( -6.16, -5.59, -5.27, -6.41, -5.74, -5.32, -6.33, -5.71, -5.33,
NA , NA, NA , -6.45, -5.67, -5.31, -6.42, -5.65, -5.32,
NA , NA, NA , NA , NA , NA , -6.32, -5.73, -5.32)
, nrow = 3, ncol = 9, byrow = TRUE )
rownames(model.two.break.tau.cval) <- c("Break 1 - 0.2", "Break 1 - 0.4", "Break 1 - 0.6")
colnames(model.two.break.tau.cval) <- c("Break 2 - 0.4 - 1%", "Break 2 - 0.4 - 5%", "Break 2 - 0.4 - 10%",
"Break 2 - 0.6 - 1%", "Break 2 - 0.6 - 5%", "Break 2 - 0.6 - 10%",
"Break 2 - 0.8 - 1%", "Break 2 - 0.8 - 5%", "Break 2 - 0.8 - 10%")
model.two.break.rho.cval <- matrix(c( -55.4 , -47.9, -44.0, -58.6, -49.9, -44.4, -57.6, -49.6, -44.6,
NA , NA, NA ,-59.3, -49.0, -44.3, -58.8, -48.7, -44.5,
NA , NA, NA , NA , NA , NA ,-57.4, -49.8, -44.4)
, nrow = 3, ncol = 9, byrow = TRUE )
rownames(model.two.break.rho.cval) <- c("Break 1 - 0.2", "Break 1 - 0.4", "Break 1 - 0.6")
colnames(model.two.break.rho.cval) <- c("Break 2 - 0.4 - 1%", "Break 2 - 0.4 - 5%", "Break 2 - 0.4 - 10%",
"Break 2 - 0.6 - 1%", "Break 2 - 0.6 - 5%", "Break 2 - 0.6 - 10%",
"Break 2 - 0.8 - 1%", "Break 2 - 0.8 - 5%", "Break 2 - 0.8 - 10%")
#Number of observations to eliminate in relation to the sample length
pnnobs <- round(pn*n)
#Define the start values
startl <- 0
myBreakStart <- startl + 1 + pnnobs
myBreakEnd <- n - pnnobs
#Calculate Dy
y.diff <- diffmatrix(y, max.diff = 1, max.lag = 1)
#Calculation
#trend for 1:n like in ur.sp
trend <- 1:n
#Define minimum gap between the two possible break dates.
#the gap is 2 in the crash case and 3 periods in the case of a break model
gap <- 2 + as.numeric(model == "break")
myBreaks <- matrix(NA, nrow = n - 2 * pnnobs, ncol = breaks)
if(breaks == 1){
myBreaks [,1] <- myBreakStart:myBreakEnd
} else if (breaks == 2){
myBreaks[, 1:breaks] <- cbind(myBreakStart:myBreakEnd,(myBreakStart:myBreakEnd)+gap)
}
#Define the variables to hold the minimum t-stat
tstat <- NA
mint <- 1000
tstat.matrix <- matrix(NA, nrow = n, ncol = n )
tstat.result <- matrix()
#Create lists to store the results
#Lists for the one break case
result.reg.coef <- list()
result.reg.resid <- list()
result.matrix <- list()
#Function to analyze the optimal lags to remove autocorrelation from the residuals
#Lag selection with general to specific procedure based on Ng,Perron (1995)
myLagSelection <- function(y.diff, S.tilde, datmat, pmax, Dummy.diff){
n <- length(y.diff)
# General to specific approach to determine the lag which removes autocorrelation from the residuals
# Ng, Perron 1995
qlag <- pmax
while (qlag >= 0) {
# Define optimal lags to include to remove autocorrelation from the residuals
# select p.value of the highest lag order and check if significant
#test.coef <- coef(summary(lm(y.diff ~ 0 + lagmatrix(S.tilde,1)[,-1] + datmat[,-1][, 1:(qlag + 1)] + Dummy.diff)))
#lm.fit implementation
test.reg.data <- na.omit(cbind(y.diff,lagmatrix(S.tilde,1)[,-1], datmat[,-1][, 1:(qlag + 1)], Dummy.diff))
test.reg.lm. <-(lm.fit(x = test.reg.data[,-1], y = test.reg.data[, 1]))
df.lm.fit <- length(test.reg.data[,1]) - test.reg.lm.$qr$rank
sigma2 <- sum((test.reg.data[,1] - test.reg.lm.$fitted.values)^2)/df.lm.fit
varbeta <- sigma2 * chol2inv(qr.R(test.reg.lm.$qr), size = ncol(test.reg.data) -2)
SE <- sqrt(diag(varbeta))
tstat <- na.omit(coef(test.reg.lm.))/SE
pvalue <- 2* pt(abs(tstat), df = df.lm.fit, lower.tail = FALSE)
# print(test.coef)
# print(paste("lm result:",qlag,test.coef[qlag + 1 , 4]))
# print(paste("lm.fit:",qlag,pvalue[qlag+1]))
# print(c("Number of qlag",qlag))
if(pvalue[qlag+1] <= 0.1){
slag <- qlag
# print("break")
break
}
qlag <- qlag - 1
slag <- qlag
}
# print(slag)
return(slag)
}
# Function to calculate the test statistic and make the code shorter, because the function can be used in both cases for
# the one break as well as the two break case
myLSfunc <- function(Dt, DTt, y.diff, est.function = c("estimation","results")){
Dt.diff <- diffmatrix(Dt, max.diff = 1, max.lag = 1)
DTt.diff <- diffmatrix(DTt, max.diff = 1, max.lag = 1)
S.tilde <- 0
S.tilde <- c(0, cumsum(lm.fit(x = na.omit(cbind(Dt.diff[,])), y=na.omit(y.diff))$residuals))
S.tilde.diff <- diffmatrix(S.tilde,max.diff = 1, max.lag = 1)
# Define optimal lags to include to remove autocorrelation
# max lag length pmax to include is based on Schwert (1989)
pmax <- min(round((12*(n/100)^0.25)),lags)
# Create matrix of lagged values of S.tilde.diff from 0 to pmax
# and check if there is autocorrelation if all these lags are included and iterate down to 1
datmat <- matrix(NA,n, pmax + 2)
datmat[ , 1] <- S.tilde.diff
# Add column of 0 to allow the easy inclusion of no lags into the test estimation
datmat[ , 2] <- rep(0, n)
if(pmax > 0){
datmat[, 3:(pmax + 2) ] <- lagmatrix(S.tilde.diff, pmax)[,-1]
colnames(datmat) <- c("S.tilde.diff", "NoLags", paste("S.tilde.diff.l",1:pmax, sep = ""))
} else if(lags == 0){
colnames(datmat) <- c("S.tilde.diff", "NoLags")
}
if(method == "Fixed"){
slag <- lags
} else if(method == "GTOS"){
slag <- NA
if(model == "crash"){
slag <- myLagSelection(y.diff, S.tilde, datmat, pmax, Dt.diff)
} else if(model == "break"){
slag <- myLagSelection(y.diff, S.tilde, datmat, pmax, DTt.diff)
}
}
S.tilde <- NA
if(model == "crash"){
S.tilde <- c(0, cumsum(lm.fit(x = na.omit(cbind(Dt.diff[,])), y=na.omit(y.diff))$residuals))
S.tilde.diff <- diffmatrix(S.tilde, max.diff = 1, max.lag = 1)
#Add lag of S.tilde.diff to control for autocorrelation in the residuals
if(est.function == "results"){
break.reg <- summary(lm(y.diff ~ 0 + lagmatrix(S.tilde, 1)[,-1] + datmat[,2:(slag+2)] + Dt.diff))
} else if (est.function == "estimation"){
#lm.fit() implementation
roll.reg.data <- na.omit(cbind(y.diff, lagmatrix(S.tilde,1)[,-1], datmat[,2:(slag+2)], Dt.diff))
roll.reg.lm. <- lm.fit(x = roll.reg.data[,-1], y = roll.reg.data[, 1])
df.lm.fit <- length(roll.reg.data[, 1]) - roll.reg.lm.$qr$rank
sigma2 <- sum((roll.reg.data[, 1] - roll.reg.lm.$fitted.values)^2)/df.lm.fit
varbeta <- sigma2*chol2inv(qr.R(roll.reg.lm.$qr), size = ncol(roll.reg.data) - 2)
SE <- sqrt(diag(varbeta))
tstat.lm.fit <- na.omit(coef(roll.reg.lm.))/SE
pvalue <- 2 * pt(abs(tstat.lm.fit),df = df.lm.fit, lower.tail = FALSE)
coef.roll.reg.lm <- cbind(na.omit(coef(roll.reg.lm.)), SE, tstat.lm.fit, pvalue)
tstat.lm.fit <- tstat.lm.fit[1]
# tstat.lm <- coef(break.reg)[1,3]
return(coef.roll.reg.lm)
}
# print(paste("lm.fit", tstat.lm.fit[1]))
# print(paste("lm", tstat.lm))
#print(roll.reg)
if(est.function == "estimation"){
return(coef.roll.reg.lm)
} else if(est.function == "results"){
return(break.reg)
}
} else if(model =="break"){
S.tilde <- c(0, cumsum(lm.fit(x = na.omit(cbind(DTt.diff[,])), y=na.omit(y.diff))$residuals))
S.tilde.diff <- diffmatrix(S.tilde, max.diff = 1, max.lag = 1)
#Add lag of S.tilde.diff to control for autocorrelation in the residuals
if(est.function == "results"){
break.reg <- summary(lm(y.diff ~ 0 + lagmatrix(S.tilde,1)[,-1] + datmat[,2:(slag+2)] + DTt.diff))
} else if (est.function == "estimation"){
#lm.fit() implementation
roll.reg.data <- na.omit(cbind(y.diff, lagmatrix(S.tilde,1)[,-1], datmat[,2:(slag+2)], DTt.diff))
roll.reg.lm. <- lm.fit(x = roll.reg.data[,-1], y = roll.reg.data[, 1])
df.lm.fit <- length(roll.reg.data[, 1]) - roll.reg.lm.$qr$rank
sigma2 <- sum((roll.reg.data[, 1] - roll.reg.lm.$fitted.values)^2)/df.lm.fit
varbeta <- sigma2*chol2inv(qr.R(roll.reg.lm.$qr), size = ncol(roll.reg.data) -2)
SE <- sqrt(diag(varbeta))
tstat.lm.fit <- na.omit(coef(roll.reg.lm.))/SE
pvalue <- 2 * pt(abs(tstat.lm.fit),df = df.lm.fit, lower.tail = FALSE)
coef.roll.reg.lm <- cbind(na.omit(coef(roll.reg.lm.)), SE, tstat.lm.fit, pvalue)
tstat.lm.fit <- tstat.lm.fit[1]
return(coef.roll.reg.lm)
# tstat.lm <- coef(break.reg)[1,3]
}
#Return Value
#print(roll.reg)
#print(paste("lm.fit", tstat.lm.fit))
#print(paste("lm", tstat.lm))
#print("break")
if(est.function == "estimation"){
return(coef.roll.reg.lm)
} else if(est.function == "results"){
return(break.reg)
}
}
# print(roll.reg)
if(est.function == "estimation"){
return(coef.roll.reg.lm)
} else if(est.function == "results"){
return(break.reg)
}
}
# Start of the actual function call
if(breaks == 1)
{
# Function to calculate the rolling t-stat
# One Break Case
for(myBreak1 in myBreaks[,1]){
#Dummies for one break case
Dt1 <- as.matrix(cbind(trend, trend >= (myBreak1 + 1)))
# Dummy with break in intercept and in trend
DTt1 <- as.matrix(cbind(Dt1, c(rep(0, myBreak1), 1:(n - myBreak1))))
colnames(Dt1) <- c("Trend","D")
colnames(DTt1) <- c("Trend","D","DTt")
#print(paste("Break1: ",myBreak1, sep = ""))
#Combine all Dummies into one big matrix to make it easier to include in the regressions
Dt <- cbind(Dt1)
DTt <- cbind(DTt1)
result.reg <- myLSfunc(Dt, DTt, y.diff, est.function = c("estimation"))
#Extract the t-statistic and if it is smaller than all previous
#t-statistics replace it and store the values of all break variables
#Extract residuals and coefficients and store them in a list
#result.matrix[[myBreak1]] <- result.reg
#result.reg.coef[[myBreak1]] <- coefficients(result.reg)
tstat <- result.reg[1,3]
tstat.result[myBreak1] <- result.reg[1,3]
#print(tstat)
if(tstat < mint){
mint <- tstat
mybestbreak1 <- myBreak1
}
}#End of first for loop
} else if(breaks == 2) {
## Two Break Case
#First for loop for the two break case
for(myBreak1 in myBreaks[,1]){
#Dummies for one break case
Dt1 <- as.matrix(cbind(trend, trend >= (myBreak1 + 1)))
# Dummy with break in intercept and in trend
DTt1 <- as.matrix(cbind(Dt1, c(rep(0, myBreak1), 1:(n - myBreak1))))
colnames(Dt1) <- c("Trend","D")
colnames(DTt1) <- c("Trend","D","DTt")
#print(paste("Break1: ",myBreak1, sep = ""))
#Second for loop for the two break case
for(myBreak2 in myBreaks[which(myBreaks[,2] < myBreakEnd & myBreaks[,2] >= myBreak1 + gap),2]){
#Dummies for two break case
Dt2 <- as.matrix(trend >= (myBreak2 + 1))
DTt2 <- as.matrix(cbind(Dt2, c(rep(0, myBreak2), 1:(n - myBreak2))))
colnames(Dt2) <- c("D2")
colnames(DTt2) <- c("D2","DTt2")
#print(paste("Break2: ",myBreak2, sep = ""))
#Combine all Dummies into one big matrix to make it easier to include in the regressions
Dt <- cbind(Dt1, Dt2)
DTt <- cbind(DTt1, DTt2)
result.reg <- myLSfunc(Dt, DTt, y.diff, est.function = c("estimation"))
#Extract the t-statistic and if it is smaller than all previous
#t-statistics replace it and store the values of all break variables
#Extract residuals and coefficients and store them in a list
#matrix to hold all the tstats
tstat.matrix[myBreak1, myBreak2] <- result.reg[1,3]
tstat <- result.reg[1,3]
#print(tstat)
if(tstat < mint){
mint <- tstat
mybestbreak1 <- myBreak1
mybestbreak2 <- myBreak2
}
}#End of second for loop
}#End of first for loop
} else if(breaks > 2){
print("Currently more than two possible structural breaks are not implemented.")
}
#Estimate regression results, based on the determined breaks and the selected lag
# to obtain all necessary information
Dt1 <- as.matrix(cbind(trend, trend >= (mybestbreak1 + 1)))
# Dummy with break in intercept and in trend
DTt1 <- as.matrix(cbind(Dt1, c(rep(0, mybestbreak1), 1:(n - mybestbreak1))))
colnames(Dt1) <- c("Trend","D")
colnames(DTt1) <- c("Trend","D","DTt")
#print(paste("Break1: ",myBreak1, sep = ""))
if(breaks == 2){
#Dummies for two break case
Dt2 <- as.matrix(trend >= (mybestbreak2 + 1))
DTt2 <- as.matrix(cbind(Dt2, c(rep(0, mybestbreak2), 1:(n - mybestbreak2))))
colnames(Dt2) <- c("D2")
colnames(DTt2) <- c("D2","DTt2")
#print(paste("Break2: ",myBreak2, sep = ""))
#Combine all Dummies into one big matrix to make it easier to include in the regressions
Dt <- cbind(Dt1, Dt2)
DTt <- cbind(DTt1, DTt2)
} else if (breaks == 1){
Dt <- Dt1
DTt <- DTt1
}
break.reg <- myLSfunc(Dt, DTt, y.diff, est.function = c("results"))
endtime <- Sys.time()
myruntime <- difftime(endtime,starttime, units = "mins")
if(print.results == "print"){
print(mint)
print(paste("First possible structural break at position:", mybestbreak1))
print(paste("The location of the first break - lambda_1:", round(mybestbreak1/n, digits = 1),", with the number of total observations:", n))
if(breaks == 2){
print(paste("Second possible structural break at position:", mybestbreak2))
print(paste("The location of the second break - lambda_2:", round(mybestbreak2/n, digits = 1),", with the number of total observations:", n))
# Output Critical Values
cat("Critical values:\n")
print(model.two.break.tau.cval)
}else if(breaks == 1){
if(model == "crash"){
cat("Critical values - Crash model:\n")
print(model.one.crash.cval)
} else if(model == "break"){
cat("Critical values - Break model:\n")
print(model.one.break.cval)
}
}
if(method == "Fixed"){
print(paste("Number of lags used:",lags))
} else if(method == "GTOS"){
print(paste("Number of lags determined by general-to-specific lag selection:"
,as.integer(substring(unlist(attr(delete.response(terms(break.reg)), "dataClasses")[3]),9))-1))
}
cat("Runtime:\n")
print(myruntime)
}
# Create complete list with all the information and not only print it
if(breaks == 2){
results.return <- list(mint, mybestbreak1, mybestbreak2, myruntime)
names(results.return) <- c("t-stat", "First break", "Second break", "Runtime")
} else if(breaks == 1){
results.return <- list(mint, mybestbreak1, myruntime)
names(results.return) <- c("t-stat", "First break", "Runtime")
}
return(list(results.return, break.reg))
}#End of ur.ls functionRunning Lee-strazicich for 1 break
# Application of the basic Lee-Strazizich test without any parallelization
# Definition of values to be used in the application
# Number of possible structural breaks
myBreaks <- 1
# Assumed break in the series, "crash" - break in intercept; "break" - break in intercept and trend
myModel <- "break"
# Number of lags to be used in fixed specification or maximum number of lags, when using the GTOS method
myLags <- 12
y <- ts_cimento
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -3.567245
## [1] "First possible structural break at position: 104"
## [1] "The location of the first break - lambda_1: 0.3 , with the number of total observations: 349"
## Critical values - Break model:
## lambda 1% 5% 10%
## [1,] 0.1 -5.11 -4.50 -4.21
## [2,] 0.2 -5.07 -4.47 -4.20
## [3,] 0.3 -5.15 -4.45 -4.18
## [4,] 0.4 -5.05 -4.50 -4.18
## [5,] 0.5 -5.11 -4.51 -4.17
## [1] "Number of lags determined by general-to-specific lag selection: 11"
## Runtime:
## Time difference of 0.02741095 minsy <- ts_calcario
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -4.867287
## [1] "First possible structural break at position: 110"
## [1] "The location of the first break - lambda_1: 0.3 , with the number of total observations: 349"
## Critical values - Break model:
## lambda 1% 5% 10%
## [1,] 0.1 -5.11 -4.50 -4.21
## [2,] 0.2 -5.07 -4.47 -4.20
## [3,] 0.3 -5.15 -4.45 -4.18
## [4,] 0.4 -5.05 -4.50 -4.18
## [5,] 0.5 -5.11 -4.51 -4.17
## [1] "Number of lags determined by general-to-specific lag selection: 12"
## Runtime:
## Time difference of 0.0114938 minsts_calcario[110]## [1] 129.4204y <- ts_diesel
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -3.876402
## [1] "First possible structural break at position: 105"
## [1] "The location of the first break - lambda_1: 0.3 , with the number of total observations: 349"
## Critical values - Break model:
## lambda 1% 5% 10%
## [1,] 0.1 -5.11 -4.50 -4.21
## [2,] 0.2 -5.07 -4.47 -4.20
## [3,] 0.3 -5.15 -4.45 -4.18
## [4,] 0.4 -5.05 -4.50 -4.18
## [5,] 0.5 -5.11 -4.51 -4.17
## [1] "Number of lags determined by general-to-specific lag selection: 7"
## Runtime:
## Time difference of 0.02056672 minsy <- ts_combustivel
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -2.903458
## [1] "First possible structural break at position: 202"
## [1] "The location of the first break - lambda_1: 0.6 , with the number of total observations: 349"
## Critical values - Break model:
## lambda 1% 5% 10%
## [1,] 0.1 -5.11 -4.50 -4.21
## [2,] 0.2 -5.07 -4.47 -4.20
## [3,] 0.3 -5.15 -4.45 -4.18
## [4,] 0.4 -5.05 -4.50 -4.18
## [5,] 0.5 -5.11 -4.51 -4.17
## [1] "Number of lags determined by general-to-specific lag selection: 9"
## Runtime:
## Time difference of 0.02870893 minsSomente a série calcário apresentou estacionariedade com uma quebra, sugerindo a quebra em Janeiro de 1989.
Running Lee-strazicich for 2 breaks
# Application of the basic Lee-Strazizich test without any parallelization
# Definition of values to be used in the application
# Number of possible structural breaks
myBreaks <- 2
# Assumed break in the series, "crash" - break in intercept; "break" - break in intercept and trend
myModel <- "break"
# Number of lags to be used in fixed specification or maximum number of lags, when using the GTOS method
myLags <- 12
y <- ts_cimento
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -4.325359
## [1] "First possible structural break at position: 84"
## [1] "The location of the first break - lambda_1: 0.2 , with the number of total observations: 349"
## [1] "Second possible structural break at position: 130"
## [1] "The location of the second break - lambda_2: 0.4 , with the number of total observations: 349"
## Critical values:
## Break 2 - 0.4 - 1% Break 2 - 0.4 - 5% Break 2 - 0.4 - 10%
## Break 1 - 0.2 -6.16 -5.59 -5.27
## Break 1 - 0.4 NA NA NA
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.6 - 1% Break 2 - 0.6 - 5% Break 2 - 0.6 - 10%
## Break 1 - 0.2 -6.41 -5.74 -5.32
## Break 1 - 0.4 -6.45 -5.67 -5.31
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.8 - 1% Break 2 - 0.8 - 5% Break 2 - 0.8 - 10%
## Break 1 - 0.2 -6.33 -5.71 -5.33
## Break 1 - 0.4 -6.42 -5.65 -5.32
## Break 1 - 0.6 -6.32 -5.73 -5.32
## [1] "Number of lags determined by general-to-specific lag selection: 8"
## Runtime:
## Time difference of 3.254675 minsy <- ts_calcario
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -5.794528
## [1] "First possible structural break at position: 102"
## [1] "The location of the first break - lambda_1: 0.3 , with the number of total observations: 349"
## [1] "Second possible structural break at position: 108"
## [1] "The location of the second break - lambda_2: 0.3 , with the number of total observations: 349"
## Critical values:
## Break 2 - 0.4 - 1% Break 2 - 0.4 - 5% Break 2 - 0.4 - 10%
## Break 1 - 0.2 -6.16 -5.59 -5.27
## Break 1 - 0.4 NA NA NA
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.6 - 1% Break 2 - 0.6 - 5% Break 2 - 0.6 - 10%
## Break 1 - 0.2 -6.41 -5.74 -5.32
## Break 1 - 0.4 -6.45 -5.67 -5.31
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.8 - 1% Break 2 - 0.8 - 5% Break 2 - 0.8 - 10%
## Break 1 - 0.2 -6.33 -5.71 -5.33
## Break 1 - 0.4 -6.42 -5.65 -5.32
## Break 1 - 0.6 -6.32 -5.73 -5.32
## [1] "Number of lags determined by general-to-specific lag selection: 11"
## Runtime:
## Time difference of 1.714561 minsy <- ts_diesel
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -4.936647
## [1] "First possible structural break at position: 42"
## [1] "The location of the first break - lambda_1: 0.1 , with the number of total observations: 349"
## [1] "Second possible structural break at position: 105"
## [1] "The location of the second break - lambda_2: 0.3 , with the number of total observations: 349"
## Critical values:
## Break 2 - 0.4 - 1% Break 2 - 0.4 - 5% Break 2 - 0.4 - 10%
## Break 1 - 0.2 -6.16 -5.59 -5.27
## Break 1 - 0.4 NA NA NA
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.6 - 1% Break 2 - 0.6 - 5% Break 2 - 0.6 - 10%
## Break 1 - 0.2 -6.41 -5.74 -5.32
## Break 1 - 0.4 -6.45 -5.67 -5.31
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.8 - 1% Break 2 - 0.8 - 5% Break 2 - 0.8 - 10%
## Break 1 - 0.2 -6.33 -5.71 -5.33
## Break 1 - 0.4 -6.42 -5.65 -5.32
## Break 1 - 0.6 -6.32 -5.73 -5.32
## [1] "Number of lags determined by general-to-specific lag selection: 8"
## Runtime:
## Time difference of 3.088442 minsy <- ts_combustivel
myLS_test <- ur.ls(y=y , model = myModel, breaks = myBreaks, lags = myLags, method = "GTOS",pn = 0.1, print.results = "print" )## [1] -4.854035
## [1] "First possible structural break at position: 107"
## [1] "The location of the first break - lambda_1: 0.3 , with the number of total observations: 349"
## [1] "Second possible structural break at position: 242"
## [1] "The location of the second break - lambda_2: 0.7 , with the number of total observations: 349"
## Critical values:
## Break 2 - 0.4 - 1% Break 2 - 0.4 - 5% Break 2 - 0.4 - 10%
## Break 1 - 0.2 -6.16 -5.59 -5.27
## Break 1 - 0.4 NA NA NA
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.6 - 1% Break 2 - 0.6 - 5% Break 2 - 0.6 - 10%
## Break 1 - 0.2 -6.41 -5.74 -5.32
## Break 1 - 0.4 -6.45 -5.67 -5.31
## Break 1 - 0.6 NA NA NA
## Break 2 - 0.8 - 1% Break 2 - 0.8 - 5% Break 2 - 0.8 - 10%
## Break 1 - 0.2 -6.33 -5.71 -5.33
## Break 1 - 0.4 -6.42 -5.65 -5.32
## Break 1 - 0.6 -6.32 -5.73 -5.32
## [1] "Number of lags determined by general-to-specific lag selection: 4"
## Runtime:
## Time difference of 3.954312 minsTestando para duas quebras os resultados se mantém, sendo apenas a série calcário considerada estacionária. Compreendendo a critica de Lee-Stravicich a Zivot-Andrews, considera-se preferivel os resultados apresentados pelo método Multiplicador de Lagrange. Desta forma compreendo que a abordagem com quebra estrutural não é muito conclusiva. Talvez seguir outra metodologia possa ser mais esclarecedor.