Instalando os pacotes
importando dados
path <- "d:/Users/thomas.daher/Desktop/cimento/dados_parte2/"
# cimento
cimento <- read_excel(paste(path, "Obras hidreletricas.xlsx", sep = ""),sheet = "cimento_51") %>%
dplyr::select(data, cimento_51_barra)
#calcário
calcario <- read_excel(paste(path, "calcario.xlsx", sep = "")) %>% dplyr::select(data,calcario)
# diesel
diesel <- read_excel(paste(path, "Obras hidreletricas.xlsx", sep = ""), sheet = "Oleo diesel") %>% dplyr::select(data,diesel)
#coque de petroleo
#coque_antigo <- read_excel("d:/Users/thomas.daher/Desktop/cimento/dados_parte2/coque de petroleo.xlsx", sheet = "dados antigos")
#coque_antigo <- coque_antigo %>% group_by(data) %>% summarise(valor = sum(valor), kg = sum(kg)) %>%mutate(valor_por_kg = valor / kg)
#coque_novo <- read_excel("d:/Users/thomas.daher/Desktop/cimento/dados_parte2/coque de petroleo.xlsx", sheet = "dados novos")
#coque_novo <- coque_novo %>% group_by(data) %>% summarise(valor = sum(`Valor FOB (US$)`), kg = sum(`Quilograma Líquido`)) %>% mutate(valor_por_kg = valor / kg)
#coque <- rbind(coque_antigo,coque_novo)
#rm(coque_antigo,coque_novo)
#producao cimento
prod_cim <- read_excel("d:/Users/thomas.daher/Desktop/cimento/dados_parte2/prod_cimento.xlsx")
prod_cim <- prod_cim %>%
pivot_longer(cols = -ano, names_to = "mes", values_to = "prod_cim") %>%
mutate(data = as.Date(paste(ano, mes, "01", sep = "-"), format = "%Y-%m-%d")) %>% dplyr::select(data, prod_cim)
# serie lubrificantes e combustiveis
combustiveis <- read_excel("d:/Users/thomas.daher/Desktop/cimento/dados_parte2/combustiveis.xlsx") %>% dplyr::select(data, combustiveis)
IPCA <- read_excel("d:/Users/thomas.daher/Desktop/cimento/dados_parte2/ipca.xlsx") %>% dplyr::select(data, IPCA_100)
Criando dataframes e deflacionando
# Criando o data.frame
df <- merge(x = merge(x = merge(x = cimento, y = diesel, by = "data", all.x = TRUE), y = calcario, by = "data", all.x = TRUE), y = combustiveis, by = "data", all.x = TRUE)
df <- df %>%
inner_join(IPCA, by = "data")
# Defina uma função para deflação
deflacionar_serie <- function(valor_nominal, ipca) {
return((valor_nominal / ipca)*100)
}
# Aplique a função de deflação para cada série no DataFrame
df <- df %>%
mutate(
cimento_deflacionado = deflacionar_serie(cimento_51_barra, IPCA_100),
diesel_deflacionado = deflacionar_serie(diesel, IPCA_100),
calcario_deflacionado = deflacionar_serie(calcario, IPCA_100),
combustiveis_deflacionado = deflacionar_serie(combustiveis, IPCA_100)
)
df <- df %>% dplyr::select(data, cimento_deflacionado, diesel_deflacionado, calcario_deflacionado, combustiveis_deflacionado, IPCA_100)
Análise de raiz unitária
# serie em nivel
adf.test(ts(df$cimento_deflacionado))
##
## Augmented Dickey-Fuller Test
##
## data: ts(df$cimento_deflacionado)
## Dickey-Fuller = -2.819, Lag order = 7, p-value = 0.2312
## alternative hypothesis: stationary
adf.test(ts(df$calcario_deflacionado))
##
## Augmented Dickey-Fuller Test
##
## data: ts(df$calcario_deflacionado)
## Dickey-Fuller = -3.6511, Lag order = 7, p-value = 0.02836
## alternative hypothesis: stationary
adf.test(ts(df$diesel_deflacionado))
##
## Augmented Dickey-Fuller Test
##
## data: ts(df$diesel_deflacionado)
## Dickey-Fuller = -1.3595, Lag order = 7, p-value = 0.847
## alternative hypothesis: stationary
adf.test(ts(df$combustiveis_deflacionado))
##
## Augmented Dickey-Fuller Test
##
## data: ts(df$combustiveis_deflacionado)
## Dickey-Fuller = -1.1255, Lag order = 7, p-value = 0.9175
## alternative hypothesis: stationary
### plotando series em nivel
ggplot(df) +
geom_line(aes(x = data, y = cimento_deflacionado, color = "cimento"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = calcario_deflacionado, color = "calcario"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = diesel_deflacionado, color = "diesel"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = combustiveis_deflacionado, color = "combustiveis"), linetype = "solid", size = 1) +
labs(title = "Séries deflacionadas de cimento, calcario, diesel e combustiveis",
x = "data",
y = "serie",
color = "Séries em nível") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# serie em primeira diferença
adf.test(ts(diff(df$cimento_deflacionado)))
## Warning in adf.test(ts(diff(df$cimento_deflacionado))): p-value smaller than
## printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts(diff(df$cimento_deflacionado))
## Dickey-Fuller = -6.4094, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ts(diff(df$calcario_deflacionado)))
## Warning in adf.test(ts(diff(df$calcario_deflacionado))): p-value smaller than
## printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts(diff(df$calcario_deflacionado))
## Dickey-Fuller = -8.6684, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ts(diff(df$diesel_deflacionado)))
## Warning in adf.test(ts(diff(df$diesel_deflacionado))): p-value smaller than
## printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts(diff(df$diesel_deflacionado))
## Dickey-Fuller = -6.5866, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ts(diff(df$combustiveis_deflacionado)))
## Warning in adf.test(ts(diff(df$combustiveis_deflacionado))): p-value smaller
## than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts(diff(df$combustiveis_deflacionado))
## Dickey-Fuller = -6.4255, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
### plotando series em primeira diferença
df$diff_cimento <- c(NA, diff(df$cimento_deflacionado))
df$diff_calcario <- c(NA, diff(df$calcario_deflacionado))
df$diff_diesel <- c(NA, diff(df$diesel_deflacionado))
df$diff_combustiveis <- c(NA, diff(df$combustiveis_deflacionado))
ggplot(df[-1, ]) + # Exclude the first row for plotting
geom_line(aes(x = data, y = diff_cimento, color = "cimento"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = diff_calcario, color = "calcario"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = diff_diesel, color = "diesel"), linetype = "solid", size = 1) +
geom_line(aes(x = data, y = diff_combustiveis, color = "combustiveis"), linetype = "solid", size = 1) +
labs(title = "séries em primeira diferença)",
x = "data",
y = "série",
color = "séries") +
theme_minimal()
# Analisando cointegração
ts_data <- ts(df[, c("diff_cimento", "diff_calcario", "diff_diesel","diff_combustiveis")], start = c(1979, 5), frequency = 12)
ts_data <- ts_data[-1, ]
#VARselect
K <- VARselect(ts_data)
K <- K$selection[1]
K
## AIC(n)
## 10
#Johansen test
johansen_result <-ca.jo(df[, c("diff_cimento", "diff_calcario", "diff_diesel", "diff_combustiveis")], type = c("eigen"), ecdet = c("none"), K = K , spec=c("longrun"), season = NULL, dumvar = NULL)
summary(johansen_result)
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.20453610 0.14254350 0.11158400 0.07985677
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 28.13 6.50 8.18 11.65
## r <= 2 | 39.99 12.91 14.90 19.19
## r <= 1 | 51.98 18.90 21.07 25.75
## r = 0 | 77.34 24.78 27.14 32.14
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## diff_cimento.l10 diff_calcario.l10 diff_diesel.l10
## diff_cimento.l10 1.0000000 1.0000000 1.00000
## diff_calcario.l10 -3.8371475 -0.7390320 57.24050
## diff_diesel.l10 1.0036844 -0.8207596 18.15873
## diff_combustiveis.l10 -0.8651031 0.9879301 -19.59666
## diff_combustiveis.l10
## diff_cimento.l10 1.00000
## diff_calcario.l10 125.60008
## diff_diesel.l10 -82.50119
## diff_combustiveis.l10 -138.72529
##
## Weights W:
## (This is the loading matrix)
##
## diff_cimento.l10 diff_calcario.l10 diff_diesel.l10
## diff_cimento.d 0.1679851 -0.55001346 -0.04518253
## diff_calcario.d 0.1920081 0.01170783 -0.02129385
## diff_diesel.d -0.9975891 0.09416912 -0.03201648
## diff_combustiveis.d 0.1818719 -0.11179055 -0.00775393
## diff_combustiveis.l10
## diff_cimento.d 0.0015813847
## diff_calcario.d 0.0008653432
## diff_diesel.d 0.0040429863
## diff_combustiveis.d 0.0049791477
MQO
years_vector <- seq(1979, 2007, by = 1)
df <- df %>% mutate(dummy_year = ifelse(year(data) >= 1987 & year(data) <= 2007, 1, 0))
formula <- "diff_cimento ~ diff_calcario + diff_diesel + diff_combustiveis + dummy_year + data"
modelo_linear <- (lm(data = df[-1, ], formula = formula))
summary(modelo_linear)
##
## Call:
## lm(formula = formula, data = df[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -76.910 -2.155 0.399 2.597 61.897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.381e-01 1.842e+00 0.455 0.650
## diff_calcario 1.861e+00 1.201e-01 15.492 <2e-16 ***
## diff_diesel 5.615e-02 3.947e-02 1.423 0.156
## diff_combustiveis 1.774e-02 6.953e-02 0.255 0.799
## dummy_year -1.462e+00 1.635e+00 -0.894 0.372
## data 3.289e-10 2.775e-09 0.119 0.906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.09 on 342 degrees of freedom
## Multiple R-squared: 0.5354, Adjusted R-squared: 0.5286
## F-statistic: 78.83 on 5 and 342 DF, p-value: < 2.2e-16
residuals_vs_time <- residuals(modelo_linear)
dates <- seq(as.Date("1980-01-01"), as.Date("2008-12-01"), by = "1 month")
plot(x = dates, y = residuals_vs_time, type = "l", xlab = "Date", ylab = "Residuals", main = "Plot dos resíduos")
MQO com prod_cimento
adf.test(ts((log(prod_cim$prod_cim))))
##
## Augmented Dickey-Fuller Test
##
## data: ts((log(prod_cim$prod_cim)))
## Dickey-Fuller = -1.7668, Lag order = 6, p-value = 0.6739
## alternative hypothesis: stationary
adf.test(ts(diff(log(prod_cim$prod_cim))))
## Warning in adf.test(ts(diff(log(prod_cim$prod_cim)))): p-value smaller than
## printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: ts(diff(log(prod_cim$prod_cim)))
## Dickey-Fuller = -10.947, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
prod_cim$diff_prod_cim <- c(NA, diff(log(prod_cim$prod_cim)))
df2 <- inner_join(x = prod_cim, y = df, by = "data")
df2 <- df2 %>% mutate(dummy_year = ifelse(year(data) >= 2003 & year(data) <= 2007, 1, 0))
formula <- "diff_cimento ~ diff_calcario + diff_diesel + diff_combustiveis + diff_prod_cim + dummy_year + data"
modelo_linear <- (lm(data = df2[-1, ], formula = formula))
summary(modelo_linear)
##
## Call:
## lm(formula = formula, data = df2[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3909 -1.1689 -0.1144 0.6944 4.5014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.176e+01 6.623e+00 -3.286 0.00165 **
## diff_calcario 1.741e+00 2.289e-01 7.608 1.58e-10 ***
## diff_diesel -1.051e-01 7.947e-02 -1.322 0.19079
## diff_combustiveis 3.837e-02 3.914e-02 0.980 0.33058
## diff_prod_cim 2.802e+00 3.151e+00 0.889 0.37731
## dummy_year 8.513e-01 7.587e-01 1.122 0.26604
## data 1.791e-08 5.495e-09 3.261 0.00178 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.801 on 64 degrees of freedom
## Multiple R-squared: 0.6341, Adjusted R-squared: 0.5998
## F-statistic: 18.48 on 6 and 64 DF, p-value: 2.497e-12
residuals_vs_time <- residuals(modelo_linear)
dates <- seq(as.Date("2003-02-01"), as.Date("2008-12-01"), by = "1 month")
plot(x = dates, y = residuals_vs_time, type = "l", xlab = "Date", ylab = "Residuals", main = "Plot dos resíduos")
Bai-Perron
ts_data <- ts(cimento$cimento_51_barra)
bp_results <- breakpoints(ts_data ~ 1)
# Extracting BIC values for different numbers of breakpoints
bic_values <- as.vector(BIC(bp_results))
# Finding the optimal number of breakpoints (minimum BIC)
optimal_breakpoints <- which.min(bic_values)
# Displaying results
summary(bp_results)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = ts_data ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 246
## m = 2 182 263
## m = 3 177 230 283
## m = 4 124 177 230 283
## m = 5 71 124 177 230 283
##
## Corresponding to breakdates:
##
## m = 1 246
## m = 2 182 263
## m = 3 177 230 283
## m = 4 124 177 230 283
## m = 5 71 124 177 230 283
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 8.605e+15 1.486e+15 3.623e+14 3.590e+14 3.589e+14 3.589e+14
## BIC 1.203e+04 1.141e+04 1.092e+04 1.093e+04 1.094e+04 1.095e+04
print(paste("Optimal number of breakpoints: ", optimal_breakpoints))
## [1] "Optimal number of breakpoints: 1"
Multivariate Bai-Perron
ts_cimento <- df[-1, ]$diff_cimento
ts_calcario <- df[-1, ]$diff_calcario
ts_combustiveis <- df[-1, ]$diff_combustiveis
# Assuming 'multivariate_data' is a data frame containing both time series
multivariate_data <- data.frame(ts_cimento, ts_calcario, ts_combustiveis)
mts_object <- ts(multivariate_data, start = c(1980, 01), frequency = 12)
bp_results_multivariate <- breakpoints(mts_object ~ 1)
# Extracting BIC values for different numbers of breakpoints
bic_values <- as.vector(BIC(bp_results_multivariate))
# Finding the optimal number of breakpoints (minimum BIC)
optimal_breakpoints <- which.min(bic_values)
# Displaying results for the multivariate time series
summary(bp_results_multivariate)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = mts_object ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 107
## m = 2 55 107
## m = 3 55 107 165
## m = 4 55 107 196 277
## m = 5 55 107 159 211 277
##
## Corresponding to breakdates:
##
## m = 1 1988(11)
## m = 2 1984(7) 1988(11)
## m = 3 1984(7) 1988(11) 1993(9)
## m = 4 1984(7) 1988(11) 1996(4) 2003(1)
## m = 5 1984(7) 1988(11) 1993(3) 1997(7) 2003(1)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 90541 89098 88024 87249 87149 87003
## BIC 2935 2941 2948 2957 2968 2979
print(paste("Optimal number of breakpoints: ", optimal_breakpoints))
## [1] "Optimal number of breakpoints: 1"
plot(multivariate_data$ts_combustiveis)
MQO colocando a dummy na quebra
years_vector <- seq(1979, 2007, by = 1)
df <- df %>%
mutate(dummy_year = case_when(
data >= as.Date("1988-11-01") ~ 1,
TRUE ~ 0
))
formula <- "diff_cimento ~ 0 + diff_calcario + diff_diesel + diff_combustiveis + dummy_year + data"
modelo_linear <- (lm(data = df[-1, ], formula = formula))
summary(modelo_linear)
##
## Call:
## lm(formula = formula, data = df[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.312 -2.202 0.165 2.927 63.683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## diff_calcario 1.841e+00 1.190e-01 15.480 <2e-16 ***
## diff_diesel 5.862e-02 3.910e-02 1.499 0.1347
## diff_combustiveis 2.304e-02 6.874e-02 0.335 0.7376
## dummy_year -5.013e+00 1.946e+00 -2.576 0.0104 *
## data 4.697e-09 1.990e-09 2.360 0.0188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.99 on 343 degrees of freedom
## Multiple R-squared: 0.5433, Adjusted R-squared: 0.5366
## F-statistic: 81.61 on 5 and 343 DF, p-value: < 2.2e-16
residuals_vs_time <- residuals(modelo_linear)
dates <- seq(as.Date("1980-01-01"), as.Date("2008-12-01"), by = "1 month")
plot(x = dates, y = residuals_vs_time, type = "l", xlab = "Date", ylab = "Residuals", main = "Plot dos resíduos")
Bai-Perron com resíduos da cointegração
#Johansen test
johansen_result <-ca.jo(df[, c("diff_cimento", "diff_calcario", "diff_diesel", "diff_combustiveis")], type = c("eigen"), ecdet = c("none"), K = K , spec=c("longrun"), season = NULL, dumvar = NULL)
summary(johansen_result)
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.20453610 0.14254350 0.11158400 0.07985677
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 28.13 6.50 8.18 11.65
## r <= 2 | 39.99 12.91 14.90 19.19
## r <= 1 | 51.98 18.90 21.07 25.75
## r = 0 | 77.34 24.78 27.14 32.14
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## diff_cimento.l10 diff_calcario.l10 diff_diesel.l10
## diff_cimento.l10 1.0000000 1.0000000 1.00000
## diff_calcario.l10 -3.8371475 -0.7390320 57.24050
## diff_diesel.l10 1.0036844 -0.8207596 18.15873
## diff_combustiveis.l10 -0.8651031 0.9879301 -19.59666
## diff_combustiveis.l10
## diff_cimento.l10 1.00000
## diff_calcario.l10 125.60008
## diff_diesel.l10 -82.50119
## diff_combustiveis.l10 -138.72529
##
## Weights W:
## (This is the loading matrix)
##
## diff_cimento.l10 diff_calcario.l10 diff_diesel.l10
## diff_cimento.d 0.1679851 -0.55001346 -0.04518253
## diff_calcario.d 0.1920081 0.01170783 -0.02129385
## diff_diesel.d -0.9975891 0.09416912 -0.03201648
## diff_combustiveis.d 0.1818719 -0.11179055 -0.00775393
## diff_combustiveis.l10
## diff_cimento.d 0.0015813847
## diff_calcario.d 0.0008653432
## diff_diesel.d 0.0040429863
## diff_combustiveis.d 0.0049791477
ts_cimento <- df[-1, ]$diff_cimento
ts_calcario <- df[-1, ]$diff_calcario
ts_combustiveis <- df[-1, ]$diff_combustiveis
# Assuming 'multivariate_data' is a data frame containing both time series
multivariate_data <- data.frame(ts_cimento, ts_calcario, ts_combustiveis)
mts_object <- ts(multivariate_data, start = c(1980, 01), frequency = 12)
K <- VARselect(multivariate_data) # optimal lag selected by AIC
K <- K$selection[1]
K
## AIC(n)
## 9
vecm_model <- VECM(
multivariate_data,
lag = K,
r = 1,
include = c("none"),
beta = NULL,
estim = c("2OLS"),
LRinclude = c("none"),
exogen = NULL
)
summary(vecm_model)
## #############
## ###Model VECM
## #############
## Full sample size: 348 End sample size: 338
## Number of variables: 3 Number of estimated slope parameters 84
## AIC 4179.037 BIC 4507.819 SSR 104600.7
## Cointegrating vector (estimated by 2OLS):
## ts_cimento ts_calcario ts_combustiveis
## r1 1 -1.88664 -0.05163684
##
##
## ECT ts_cimento -1
## Equation ts_cimento -0.4362(0.2135)* -0.1210(0.2018)
## Equation ts_calcario 0.2131(0.0916)* -0.2087(0.0866)*
## Equation ts_combustiveis 0.2337(0.1752) -0.2160(0.1656)
## ts_calcario -1 ts_combustiveis -1
## Equation ts_cimento -1.7304(0.3916)*** -0.1623(0.0801)*
## Equation ts_calcario -0.6844(0.1680)*** -0.0235(0.0344)
## Equation ts_combustiveis 0.2034(0.3213) -1.0168(0.0657)***
## ts_cimento -2 ts_calcario -2
## Equation ts_cimento 0.0723(0.1950) -2.0738(0.3890)***
## Equation ts_calcario -0.0976(0.0837) -0.7786(0.1669)***
## Equation ts_combustiveis -0.1296(0.1600) 0.0841(0.3191)
## ts_combustiveis -2 ts_cimento -3
## Equation ts_cimento -0.4126(0.1115)*** -0.1523(0.1893)
## Equation ts_calcario -0.0892(0.0478). -0.1316(0.0812)
## Equation ts_combustiveis -0.8758(0.0914)*** -0.1580(0.1553)
## ts_calcario -3 ts_combustiveis -3
## Equation ts_cimento -1.9649(0.3914)*** -0.5308(0.1260)***
## Equation ts_calcario -0.7271(0.1679)*** -0.1871(0.0541)***
## Equation ts_combustiveis 0.1950(0.3211) -0.8336(0.1034)***
## ts_cimento -4 ts_calcario -4
## Equation ts_cimento 0.1559(0.1793) -1.6980(0.3873)***
## Equation ts_calcario 0.0356(0.0769) -0.6981(0.1662)***
## Equation ts_combustiveis 0.0426(0.1471) 0.1768(0.3178)
## ts_combustiveis -4 ts_cimento -5
## Equation ts_cimento -0.4858(0.1339)*** 0.1110(0.1606)
## Equation ts_calcario -0.2025(0.0574)*** -0.0259(0.0689)
## Equation ts_combustiveis -0.7863(0.1098)*** -0.1155(0.1318)
## ts_calcario -5 ts_combustiveis -5
## Equation ts_cimento -1.4683(0.3653)*** -0.3625(0.1393)**
## Equation ts_calcario -0.5639(0.1567)*** -0.1486(0.0598)*
## Equation ts_combustiveis 0.3525(0.2997) -0.6177(0.1143)***
## ts_cimento -6 ts_calcario -6
## Equation ts_cimento 0.3176(0.1407)* -1.2069(0.3252)***
## Equation ts_calcario 0.0069(0.0604) -0.4413(0.1395)**
## Equation ts_combustiveis -0.0443(0.1155) 0.5032(0.2668).
## ts_combustiveis -6 ts_cimento -7
## Equation ts_cimento -0.2806(0.1373)* 0.4727(0.1173)***
## Equation ts_calcario -0.0814(0.0589) 0.0776(0.0503)
## Equation ts_combustiveis -0.5200(0.1126)*** 0.0013(0.0963)
## ts_calcario -7 ts_combustiveis -7
## Equation ts_cimento -1.0012(0.2881)*** -0.0783(0.1314)
## Equation ts_calcario -0.3873(0.1236)** -0.0349(0.0564)
## Equation ts_combustiveis 0.4940(0.2364)* -0.3830(0.1078)***
## ts_cimento -8 ts_calcario -8
## Equation ts_cimento 0.2283(0.1035)* -0.8117(0.2528)**
## Equation ts_calcario 0.0145(0.0444) -0.3497(0.1085)**
## Equation ts_combustiveis 0.0230(0.0849) 0.1084(0.2074)
## ts_combustiveis -8 ts_cimento -9
## Equation ts_cimento 0.0539(0.1165) 0.1622(0.0920).
## Equation ts_calcario 0.0002(0.0500) 0.0288(0.0395)
## Equation ts_combustiveis -0.1784(0.0956). -0.0586(0.0755)
## ts_calcario -9 ts_combustiveis -9
## Equation ts_cimento -0.6451(0.2043)** 0.0815(0.0808)
## Equation ts_calcario -0.2134(0.0877)* 0.0110(0.0347)
## Equation ts_combustiveis 0.2176(0.1677) -0.1264(0.0663).
residuals_vecm <- residuals(vecm_model)
residuals_vecm <- ts(residuals_vecm, start = c(1980, 01), frequency = 12)
bp_results_multivariate <- breakpoints(residuals_vecm ~ 1)
# Extracting BIC values for different numbers of breakpoints
bic_values <- as.vector(BIC(bp_results_multivariate))
# Finding the optimal number of breakpoints (minimum BIC)
optimal_breakpoints <- which.min(bic_values)
# Displaying results for the multivariate time series
summary(bp_results_multivariate)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = residuals_vecm ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 97
## m = 2 59 109
## m = 3 59 109 175
## m = 4 59 109 177 268
## m = 5 59 109 164 218 268
##
## Corresponding to breakdates:
##
## m = 1 1988(1)
## m = 2 1984(11) 1989(1)
## m = 3 1984(11) 1989(1) 1994(7)
## m = 4 1984(11) 1989(1) 1994(9) 2002(4)
## m = 5 1984(11) 1989(1) 1993(8) 1998(2) 2002(4)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 56322 56217 55881 55724 55627 55633
## BIC 2700 2711 2721 2731 2742 2754
print(paste("Optimal number of breakpoints: ", optimal_breakpoints))
## [1] "Optimal number of breakpoints: 1"
Testes deixados de lado
Chow test
# plotando cada serie
layout(matrix(1:3, nrow = 3, ncol = 1))
plot(df$diff_cimento, type = "l", xlab = "Time", ylab = "cimento")
plot(df$diff_calcario, type = "l", xlab = "Time", ylab = "calcario")
plot(df$diff_combustiveis, type = "l", xlab = "Time", ylab = "combustiveis")
plot(x = df$diff_cimento, y = df$diff_calcario, xlab = "cimento", ylab = "calcario")
sctest(diff_cimento ~ diff_calcario , type = "Chow", data = df)
##
## Chow test
##
## data: diff_cimento ~ diff_calcario
## F = 4.4954, p-value = 0.01182
plot(x = df$diff_cimento, y = df$diff_combustiveis, xlab = "cimento", ylab = "combustiveis")
sctest(diff_cimento ~ diff_combustiveis, type = "Chow", data = df)
##
## Chow test
##
## data: diff_cimento ~ diff_combustiveis
## F = 6.236, p-value = 0.002186
sctest(diff_cimento ~ diff_combustiveis + diff_calcario, type = "Chow", data = df)
##
## Chow test
##
## data: diff_cimento ~ diff_combustiveis + diff_calcario
## F = 2.9108, p-value = 0.03456
#Ploberger and Kr¨amer (1992) suggested to base a structural change test on cumulative sums of the common OLS residuals. Thus, the OLS-CUSUM type empirical fluctuation process is defined by:
efp test
efp_test <- efp(diff_cimento ~ diff_calcario + diff_combustiveis + dummy_year, type = "OLS-CUSUM", data = df, lrvar = FALSE)
plot(efp_test)
# Check for significant breakpoints
if (!is.null(efp_test$breakpoints)) {
# Display the breakpoints
print(efp_test$breakpoints)
# Display the CUSUM test statistics
if (!is.null(efp_test$cusum)) {
print(efp_test$cusum)
} else {
cat("CUSUM test statistics not available.\n")
}
} else {
cat("No significant breakpoints detected.\n")
}
## No significant breakpoints detected.