#========= VAR PADRÃO

dados_var <- tibble(
  log_despesa = log(as.numeric(despesa_real_tri_sa)),
  log_pib     = log(as.numeric(pib_real_tri_sa)),
  log_receita = log(as.numeric(receita_real_tri_sa))
)

#Seleção do número de defasagens
lagselect <- VARselect(dados_var, lag.max = 8, type = "const")
lagselect
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

$criteria
                         1                   2                   3
AIC(n) -18.507288951581771 -18.434069981855131 -18.437546470750554
HQ(n)  -18.374530442417939 -18.201742590818426 -18.105650197840973
SC(n)  -18.178360006749287 -17.858444328398285 -17.615224108669345
FPE(n)   0.000000009171872   0.000000009875727   0.000000009858225
                         4                  5                  6
AIC(n) -18.484916431208049 -18.34479451817112 -18.25451908393994
HQ(n)  -18.053451276425598 -17.81376048151579 -17.62391616541174
SC(n)  -17.415897360502480 -17.02907873884119 -16.69210659598564
FPE(n)   0.000000009431553   0.00000001090462   0.00000001202309
                        7                 8
AIC(n) -18.13277085431729 -18.0442624520824
HQ(n)  -17.40259905391622 -17.2145217698084
SC(n)  -16.32366165773863 -15.9884565468793
FPE(n)   0.00000001371955   0.0000000151957
#Estimação do VAR
var_model <- VAR(dados_var, p = 1, type = "const")
summary(var_model)

VAR Estimation Results:
========================= 
Endogenous variables: log_despesa, log_pib, log_receita 
Deterministic variables: const 
Sample size: 99 
Log Likelihood: 509.768 
Roots of the characteristic polynomial:
0.9889 0.4696 0.2294
Call:
VAR(y = dados_var, p = 1, type = "const")


Estimation results for equation log_despesa: 
============================================ 
log_despesa = log_despesa.l1 + log_pib.l1 + log_receita.l1 + const 

               Estimate Std. Error t value    Pr(>|t|)    
log_despesa.l1  0.47235    0.08408   5.618 0.000000192 ***
log_pib.l1      1.03351    0.19294   5.357 0.000000590 ***
log_receita.l1 -0.35136    0.11422  -3.076     0.00274 ** 
const          -3.74785    0.77157  -4.857 0.000004676 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.07609 on 95 degrees of freedom
Multiple R-Squared: 0.9558, Adjusted R-squared: 0.9544 
F-statistic: 684.1 on 3 and 95 DF,  p-value: < 0.00000000000000022 


Estimation results for equation log_pib: 
======================================== 
log_pib = log_despesa.l1 + log_pib.l1 + log_receita.l1 + const 

               Estimate Std. Error t value             Pr(>|t|)    
log_despesa.l1  0.06859    0.02081   3.297             0.001377 ** 
log_pib.l1      0.86661    0.04774  18.151 < 0.0000000000000002 ***
log_receita.l1  0.03155    0.02826   1.116             0.267132    
const           0.66405    0.19093   3.478             0.000764 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.01883 on 95 degrees of freedom
Multiple R-Squared: 0.9949, Adjusted R-squared: 0.9948 
F-statistic:  6237 on 3 and 95 DF,  p-value: < 0.00000000000000022 


Estimation results for equation log_receita: 
============================================ 
log_receita = log_despesa.l1 + log_pib.l1 + log_receita.l1 + const 

               Estimate Std. Error t value Pr(>|t|)   
log_despesa.l1  0.03366    0.08974   0.375  0.70840   
log_pib.l1      0.58765    0.20592   2.854  0.00530 **
log_receita.l1  0.34898    0.12190   2.863  0.00517 **
const          -0.61773    0.82346  -0.750  0.45501   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.08121 on 95 degrees of freedom
Multiple R-Squared: 0.9135, Adjusted R-squared: 0.9107 
F-statistic: 334.3 on 3 and 95 DF,  p-value: < 0.00000000000000022 



Covariance matrix of residuals:
            log_despesa    log_pib log_receita
log_despesa   0.0057899 -0.0002934  -0.0016205
log_pib      -0.0002934  0.0003545   0.0009564
log_receita  -0.0016205  0.0009564   0.0065949

Correlation matrix of residuals:
            log_despesa log_pib log_receita
log_despesa      1.0000 -0.2048     -0.2623
log_pib         -0.2048  1.0000      0.6254
log_receita     -0.2623  0.6254      1.0000
#Funções de Impulso-Resposta (IRF)

#Ordem de Granger

irf_var_padrão_pib <- irf(var_model,
                          impulse = c("log_pib", "log_despesa", "log_receita"),
                          response = c("log_pib"),
                          n.ahead = 36,     
                          boot = TRUE)     
plot(irf_var_padrão_pib)

#Testes de robustez / Diagnóstico

#Teste de autocorrelação dos resíduos (LM test)
serial.test(var_model, lags.pt = 16, type = "PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var_model
Chi-squared = 121.74, df = 135, p-value = 0.7864
#Teste de normalidade dos resíduos
normality.test(var_model)
$JB

    JB-Test (multivariate)

data:  Residuals of VAR object var_model
Chi-squared = 2546.1, df = 6, p-value < 0.00000000000000022


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object var_model
Chi-squared = 252.81, df = 3, p-value < 0.00000000000000022


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object var_model
Chi-squared = 2293.3, df = 3, p-value < 0.00000000000000022
## Teste de heterocedasticidade (ARCH LM test)
arch.test(var_model, lags.multi = 12)

    ARCH (multivariate)

data:  Residuals of VAR object var_model
Chi-squared = 449.31, df = 432, p-value = 0.273
#Estabilidade do modelo (roots dentro do círculo unitário)
roots(var_model) #estabilidade
[1] 0.9889224 0.4696088 0.2294131
#========= VAR EM 1° DIFERENÇA

dados_var_diff <- tibble(
  d_log_despesa = diff(log(as.numeric(despesa_real_tri_sa))),
  d_log_pib     = diff(log(as.numeric(pib_real_tri_sa))),
  d_log_receita = diff(log(as.numeric(receita_real_tri_sa)))
)

#Seleção do número de defasagens
lagselect <- VARselect(dados_var_diff, lag.max = 8, type = "const")
lagselect
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      3      1      3 

$criteria
                       1                  2                  3
AIC(n) -17.8988895880874 -18.08261220302435 -18.24033165289189
HQ(n)  -17.7653104109537 -17.84884864304039 -17.90638371005766
SC(n)  -17.5677872355797 -17.50318308613584 -17.41257577162260
FPE(n)   0.0000000168533   0.00000001403518   0.00000001200841
                        4                 5                  6
AIC(n) -18.13110412659157 -18.0582788321250 -17.97429299993324
HQ(n)  -17.69697180090708 -17.5239621235902 -17.33979190854821
SC(n)  -17.05502148094149 -16.7338694220941 -16.40155682552158
FPE(n)   0.00000001343771   0.0000000145278   0.00000001592147
                        7                  8
AIC(n) -17.86902035500920 -17.80495499525326
HQ(n)  -17.13433488077391 -16.97008513816771
SC(n)  -16.04795741621676 -15.73556529208004
FPE(n)   0.00000001787745   0.00000001933207
#Estimação do VAR em diferenças
var_diff <- VAR(dados_var_diff, p = 3, type = "const")
summary(var_diff)

VAR Estimation Results:
========================= 
Endogenous variables: d_log_despesa, d_log_pib, d_log_receita 
Deterministic variables: const 
Sample size: 96 
Log Likelihood: 500.544 
Roots of the characteristic polynomial:
0.789 0.7527 0.7527 0.6704 0.6704 0.4681 0.4681 0.2145 0.2145
Call:
VAR(y = dados_var_diff, p = 3, type = "const")


Estimation results for equation d_log_despesa: 
============================================== 
d_log_despesa = d_log_despesa.l1 + d_log_pib.l1 + d_log_receita.l1 + d_log_despesa.l2 + d_log_pib.l2 + d_log_receita.l2 + d_log_despesa.l3 + d_log_pib.l3 + d_log_receita.l3 + const 

                 Estimate Std. Error t value Pr(>|t|)    
d_log_despesa.l1 -0.38167    0.11086  -3.443 0.000892 ***
d_log_pib.l1     -0.12232    0.54366  -0.225 0.822517    
d_log_receita.l1 -0.25980    0.12917  -2.011 0.047429 *  
d_log_despesa.l2 -0.15511    0.11545  -1.343 0.182657    
d_log_pib.l2     -0.18091    0.52661  -0.344 0.732032    
d_log_receita.l2  0.09702    0.14498   0.669 0.505135    
d_log_despesa.l3 -0.07656    0.10392  -0.737 0.463328    
d_log_pib.l3      0.36893    0.52503   0.703 0.484145    
d_log_receita.l3  0.17921    0.13405   1.337 0.184780    
const             0.01703    0.01089   1.563 0.121706    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.08028 on 86 degrees of freedom
Multiple R-Squared: 0.2892, Adjusted R-squared: 0.2148 
F-statistic: 3.887 on 9 and 86 DF,  p-value: 0.0003604 


Estimation results for equation d_log_pib: 
========================================== 
d_log_pib = d_log_despesa.l1 + d_log_pib.l1 + d_log_receita.l1 + d_log_despesa.l2 + d_log_pib.l2 + d_log_receita.l2 + d_log_despesa.l3 + d_log_pib.l3 + d_log_receita.l3 + const 

                  Estimate Std. Error t value Pr(>|t|)   
d_log_despesa.l1  0.050754   0.024876   2.040  0.04439 * 
d_log_pib.l1      0.160074   0.121989   1.312  0.19294   
d_log_receita.l1 -0.011885   0.028984  -0.410  0.68278   
d_log_despesa.l2  0.034949   0.025906   1.349  0.18085   
d_log_pib.l2      0.372693   0.118161   3.154  0.00222 **
d_log_receita.l2 -0.075273   0.032530  -2.314  0.02306 * 
d_log_despesa.l3  0.048134   0.023318   2.064  0.04201 * 
d_log_pib.l3     -0.272540   0.117807  -2.313  0.02308 * 
d_log_receita.l3 -0.025130   0.030078  -0.835  0.40576   
const             0.005952   0.002445   2.435  0.01698 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.01801 on 86 degrees of freedom
Multiple R-Squared: 0.2727, Adjusted R-squared: 0.1965 
F-statistic: 3.582 on 9 and 86 DF,  p-value: 0.0008013 


Estimation results for equation d_log_receita: 
============================================== 
d_log_receita = d_log_despesa.l1 + d_log_pib.l1 + d_log_receita.l1 + d_log_despesa.l2 + d_log_pib.l2 + d_log_receita.l2 + d_log_despesa.l3 + d_log_pib.l3 + d_log_receita.l3 + const 

                  Estimate Std. Error t value   Pr(>|t|)    
d_log_despesa.l1  0.006143   0.115095   0.053     0.9576    
d_log_pib.l1      0.994616   0.564418   1.762     0.0816 .  
d_log_receita.l1 -0.653731   0.134101  -4.875 0.00000493 ***
d_log_despesa.l2 -0.008694   0.119860  -0.073     0.9423    
d_log_pib.l2      1.389567   0.546709   2.542     0.0128 *  
d_log_receita.l2 -0.712865   0.150511  -4.736 0.00000851 ***
d_log_despesa.l3 -0.017875   0.107890  -0.166     0.8688    
d_log_pib.l3     -0.079625   0.545071  -0.146     0.8842    
d_log_receita.l3 -0.252399   0.139163  -1.814     0.0732 .  
const             0.005817   0.011311   0.514     0.6084    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.08335 on 86 degrees of freedom
Multiple R-Squared: 0.3182, Adjusted R-squared: 0.2468 
F-statistic: 4.459 on 9 and 86 DF,  p-value: 0.00008169 



Covariance matrix of residuals:
              d_log_despesa  d_log_pib d_log_receita
d_log_despesa     0.0064457 -0.0003081    -0.0019063
d_log_pib        -0.0003081  0.0003245     0.0009252
d_log_receita    -0.0019063  0.0009252     0.0069471

Correlation matrix of residuals:
              d_log_despesa d_log_pib d_log_receita
d_log_despesa        1.0000   -0.2130       -0.2849
d_log_pib           -0.2130    1.0000        0.6162
d_log_receita       -0.2849    0.6162        1.0000
#Impulso-resposta

irf_diff_granger_pib <- irf(var_diff,
                            impulse  = c("d_log_pib", "d_log_despesa", "d_log_receita"),
                            response = c("d_log_pib"),
                            n.ahead = 36,
                            boot = TRUE)
plot(irf_diff_granger_pib)

#Testes de robustez

#Autocorrelação dos resíduos
serial.test(var_diff, lags.pt = 16, type = "PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var_diff
Chi-squared = 98.014, df = 117, p-value = 0.8982
#Normalidade dos resíduos
normality.test(var_diff)
$JB

    JB-Test (multivariate)

data:  Residuals of VAR object var_diff
Chi-squared = 751.54, df = 6, p-value < 0.00000000000000022


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object var_diff
Chi-squared = 105.76, df = 3, p-value < 0.00000000000000022


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object var_diff
Chi-squared = 645.78, df = 3, p-value < 0.00000000000000022
#Heterocedasticidade (ARCH)
arch.test(var_diff, lags.multi = 12)

    ARCH (multivariate)

data:  Residuals of VAR object var_diff
Chi-squared = 434.33, df = 432, p-value = 0.4595
## Estabilidade do VAR
roots(var_diff) 
[1] 0.7889864 0.7526535 0.7526535 0.6703958 0.6703958 0.4681112 0.4681112
[8] 0.2144888 0.2144888
#========= VECM

#========= TESTE DE COINTEGRAÇÃO (Johansen)

# Determinar número de defasagens (p) para o VECM
lagselect <- VARselect(dados_var, lag.max = 8, type = "const")
lagselect
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      1 

$criteria
                         1                   2                   3
AIC(n) -18.507288951581771 -18.434069981855131 -18.437546470750554
HQ(n)  -18.374530442417939 -18.201742590818426 -18.105650197840973
SC(n)  -18.178360006749287 -17.858444328398285 -17.615224108669345
FPE(n)   0.000000009171872   0.000000009875727   0.000000009858225
                         4                  5                  6
AIC(n) -18.484916431208049 -18.34479451817112 -18.25451908393994
HQ(n)  -18.053451276425598 -17.81376048151579 -17.62391616541174
SC(n)  -17.415897360502480 -17.02907873884119 -16.69210659598564
FPE(n)   0.000000009431553   0.00000001090462   0.00000001202309
                        7                 8
AIC(n) -18.13277085431729 -18.0442624520824
HQ(n)  -17.40259905391622 -17.2145217698084
SC(n)  -16.32366165773863 -15.9884565468793
FPE(n)   0.00000001371955   0.0000000151957
# Rodar teste de Johansen
vec_test <- ca.jo(dados_var, type = "trace", K = 2, ecdet = "const", spec = "transitory")
summary(vec_test)  # número de relações de cointegração

###################### 
# Johansen-Procedure # 
###################### 

Test type: trace statistic , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 0.328575247205597387179 0.214399621246646548967 0.173171171238841348217
[4] 0.000000000000003996803

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 2 | 18.64  7.52  9.24 12.97
r <= 1 | 42.28 17.85 19.96 24.60
r = 0  | 81.32 32.00 34.91 41.07

Eigenvectors, normalised to first column:
(These are the cointegration relations)

               log_despesa.l1  log_pib.l1 log_receita.l1  constant
log_despesa.l1       1.000000  1.00000000       1.000000   1.00000
log_pib.l1          -3.445401 -1.29087188       1.144834 -10.30656
log_receita.l1       2.084027  0.04471097      -2.926712   2.43120
constant            10.605893  5.32731660       8.364458 104.99142

Weights W:
(This is the loading matrix)

              log_despesa.l1  log_pib.l1 log_receita.l1
log_despesa.d    -0.03924814 -0.37630156   -0.004812536
log_pib.d         0.04115290  0.01145064    0.020458894
log_receita.d    -0.23246106  0.11014225    0.079202754
                              constant
log_despesa.d  0.000000000000037723893
log_pib.d     -0.000000000000006082256
log_receita.d  0.000000000000034495979
# Rodar teste de Johansen
vec_test <- ca.jo(dados_var, type = "eigen", K = 2, ecdet = "const", spec = "transitory")
summary(vec_test)

###################### 
# Johansen-Procedure # 
###################### 

Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration 

Eigenvalues (lambda):
[1] 0.328575247205597387179 0.214399621246646548967 0.173171171238841348217
[4] 0.000000000000003996803

Values of teststatistic and critical values of test:

          test 10pct  5pct  1pct
r <= 2 | 18.64  7.52  9.24 12.97
r <= 1 | 23.65 13.75 15.67 20.20
r = 0  | 39.04 19.77 22.00 26.81

Eigenvectors, normalised to first column:
(These are the cointegration relations)

               log_despesa.l1  log_pib.l1 log_receita.l1  constant
log_despesa.l1       1.000000  1.00000000       1.000000   1.00000
log_pib.l1          -3.445401 -1.29087188       1.144834 -10.30656
log_receita.l1       2.084027  0.04471097      -2.926712   2.43120
constant            10.605893  5.32731660       8.364458 104.99142

Weights W:
(This is the loading matrix)

              log_despesa.l1  log_pib.l1 log_receita.l1
log_despesa.d    -0.03924814 -0.37630156   -0.004812536
log_pib.d         0.04115290  0.01145064    0.020458894
log_receita.d    -0.23246106  0.11014225    0.079202754
                              constant
log_despesa.d  0.000000000000037723893
log_pib.d     -0.000000000000006082256
log_receita.d  0.000000000000034495979
# Estimar VECM
vec_mod <- cajorls(vec_test, r = 2)
summary(vec_mod$rlm)
Response log_despesa.d :

Call:
lm(formula = log_despesa.d ~ ect1 + ect2 + log_despesa.dl1 + 
    log_pib.dl1 + log_receita.dl1 - 1, data = data.mat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.15240 -0.02510 -0.00669  0.01163  0.46098 

Coefficients:
                Estimate Std. Error t value  Pr(>|t|)    
ect1            -0.41555    0.10028  -4.144 0.0000752 ***
ect2             0.62098    0.24856   2.498    0.0142 *  
log_despesa.dl1 -0.10275    0.09685  -1.061    0.2914    
log_pib.dl1     -0.21170    0.44732  -0.473    0.6371    
log_receita.dl1 -0.27439    0.11486  -2.389    0.0189 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07367 on 93 degrees of freedom
Multiple R-squared:  0.3633,    Adjusted R-squared:  0.3291 
F-statistic: 10.61 on 5 and 93 DF,  p-value: 0.00000004381


Response log_pib.d :

Call:
lm(formula = log_pib.d ~ ect1 + ect2 + log_despesa.dl1 + log_pib.dl1 + 
    log_receita.dl1 - 1, data = data.mat)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.088441 -0.005552  0.006162  0.014484  0.076099 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
ect1             0.05260    0.02800   1.879   0.0634 .
ect2            -0.15657    0.06940  -2.256   0.0264 *
log_despesa.dl1  0.02960    0.02704   1.094   0.2766  
log_pib.dl1      0.23444    0.12489   1.877   0.0636 .
log_receita.dl1 -0.03426    0.03207  -1.068   0.2882  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02057 on 93 degrees of freedom
Multiple R-squared:  0.1419,    Adjusted R-squared:  0.09578 
F-statistic: 3.076 on 5 and 93 DF,  p-value: 0.01295


Response log_receita.d :

Call:
lm(formula = log_receita.d ~ ect1 + ect2 + log_despesa.dl1 + 
    log_pib.dl1 + log_receita.dl1 - 1, data = data.mat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.47344 -0.01317  0.01951  0.05334  0.39693 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
ect1             -0.1223     0.1179  -1.038   0.3022  
ect2              0.6587     0.2922   2.254   0.0265 *
log_despesa.dl1   0.2912     0.1139   2.558   0.0121 *
log_pib.dl1       0.8589     0.5259   1.633   0.1058  
log_receita.dl1  -0.1145     0.1350  -0.848   0.3985  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0866 on 93 degrees of freedom
Multiple R-squared:  0.2134,    Adjusted R-squared:  0.1712 
F-statistic: 5.047 on 5 and 93 DF,  p-value: 0.0003888
# Transformar o resultado do ca.jo em um objeto VAR/VEC compatível
vec2var_mod <- vec2var(vec_test, r = 2)

#Testes de robustez / diagnóstico

# Autocorrelação dos resíduos (Portmanteau)
serial.test(vec2var_mod, lags.pt = 12, type = "PT.asymptotic")

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object vec2var_mod
Chi-squared = 91.227, df = 93, p-value = 0.5326
# Normalidade dos resíduos
normality.test(vec2var_mod)
$JB

    JB-Test (multivariate)

data:  Residuals of VAR object vec2var_mod
Chi-squared = 2048.2, df = 6, p-value < 0.00000000000000022


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object vec2var_mod
Chi-squared = 233.16, df = 3, p-value < 0.00000000000000022


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object vec2var_mod
Chi-squared = 1815.1, df = 3, p-value < 0.00000000000000022
# Heterocedasticidade (ARCH)
arch.test(vec2var_mod, lags.multi = 12)

    ARCH (multivariate)

data:  Residuals of VAR object vec2var_mod
Chi-squared = 458.06, df = 432, p-value = 0.1862
#Impulso-resposta

irf_despesa_pib <- irf(vec2var_mod, 
                       impulse = "log_despesa", 
                       response = "log_pib", 
                       n.ahead = 36, 
                       ci = 0.95, 
                       boot = TRUE)

plot(irf_despesa_pib)

# Calcule o impulso-resposta de um choque na receita sobre o pib
irf_receita_pib <- irf(vec2var_mod, 
                       impulse = "log_receita", 
                       response = "log_pib", 
                       n.ahead = 36, 
                       ci = 0.95, 
                       boot = TRUE)

plot(irf_receita_pib)