####################################################################################
###AULA 17/07/23
###################################################################################
###################################################################################
#limpando a memoria
rm(list=ls(all=TRUE)) 
####################################################################################
####################################################################################
#Pacotes necessários
require(tseries)
require(urca)
require(vars)
require(gap)
require(fUnitRoots)
require(MSBVAR)
require(dynlm) 
require(stats) 
library(forecast) 
require(strucchange)
#####################################################################
#Explorando o pacote vars
#####################################################################
##carregando a base de dados
par(mfrow=c(1,1))
data(Canada)
?Canada
plot(Canada)
summary(Canada)
Canada
emprego = (Canada[,"e"])
emprego
length(emprego)
Canada= data.frame(Canada)
names(Canada)
#################################################################################
#series do arquivo Canadá
#e = emprego
#prod = produtividade do trabalho
#rw = salário real
#U = taxa de desemprego

##############################################################################
##Teste de raíz unitária ADF utilizando o pacote urca 
##############################################################################
## dados emprego

adfet <- summary(ur.df(Canada[,"e"], type="trend", lags=2))
adfet

adfec <- summary(ur.df(Canada[,"e"], type="drift", selectlags=c("AIC")))
adfec

adfen <- summary(ur.df(Canada[,"e"], type="none", selectlags=c("AIC")))
adfen

adfedif<- summary(ur.df(diff(Canada[,"e"]), type="none", selectlags=c("AIC")))
adfedif

################################################################################
################################################################################
## dados salario real

adfrwt <- summary(ur.df(Canada[,"rw"], type="trend", lags=2))
adfrwt

adfrwc <- summary(ur.df(Canada[,"rw"], type="drift", selectlags=c("AIC")))
adfrwc

adfrwn <- summary(ur.df(Canada[,"rw"], type="none", selectlags=c("AIC")))
adfrwn

adfrwdif<- summary(ur.df(diff(Canada[,"rw"]), type="none", selectlags=c("AIC")))
adfrwdif

################################################################################
## dados produtividade

adfprodt <- summary(ur.df(Canada[,"prod"], type="trend", lags=2))
adfprodt

adfprodc <- summary(ur.df(Canada[,"prod"], type="drift", selectlags=c("AIC")))
adfprodc

adfprodn <- summary(ur.df(Canada[,"prod"], type="none", selectlags=c("AIC")))
adfprodn

adfproddif<- summary(ur.df(diff(Canada[,"prod"]), type="none", selectlags=c("AIC")))
adfproddif
################################################################################
## dados taxa de desmprego

adfUt <- summary(ur.df(Canada[,"U"], type="trend", lags=2))
adfUt

adfUc <- summary(ur.df(Canada[,"U"], type="drift", selectlags=c("AIC")))
adfUc

adfUn <- summary(ur.df(Canada[,"U"], type="none", selectlags=c("AIC")))
adfUn

adfUdif<- summary(ur.df(diff(Canada[,"U"]), type="none", selectlags=c("AIC")))
adfUdif

##############################################################################
# teste DF-GLS 
#SERIE emprego - e

DFglset= ur.ers(Canada[,"e"], type= "DF-GLS", model="trend", lag.max=4)
summary(DFglset)

DFglsdet= ur.ers(diff(Canada[,"e"]), type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsdet)

DFglsec= ur.ers(Canada[,"e"], type= "DF-GLS", model="const", lag.max=4)
summary(DFglsec)


DFglsdec= ur.ers(diff(Canada[,"e"]), type= "DF-GLS", model="const", lag.max=4)
summary(DFglsdec)

#########################################################################################
# teste DF-GLS 
#SERIE salário real - rw
DFglsrwt= ur.ers(Canada[,"rw"], type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsrwt)

DFglsdrwt= ur.ers(diff(Canada[,"rw"]), type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsdrwt)

DFglsdrwc= ur.ers(Canada[,"rw"], type= "DF-GLS", model="const", lag.max=4)
summary(DFglsrwc)

DFglsdrwc= ur.ers(diff(Canada[,"rw"]), type= "DF-GLS", model="const", lag.max=4)
summary(DFglsdrwc)

#########################################################################################
# teste DF-GLS 
#SERIE produtividade - prod
DFglsprodt= ur.ers(Canada[,"prod"], type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsprodt)

DFglsdprodt= ur.ers(diff(Canada[,"prod"]), type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsdprodt)

DFglsdprodc= ur.ers(Canada[,"prod"], type= "DF-GLS", model="const", lag.max=4)
summary(DFglsprodc)

DFglsdprodc= ur.ers(diff(Canada[,"prod"]), type= "DF-GLS", model="const", lag.max=4)
summary(DFglsdprodc)

#########################################################################################
# teste DF-GLS 
#SERIE desemprego - U
DFglsUt= ur.ers(Canada[,"U"], type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsUt)

DFglsdUt= ur.ers(diff(Canada[,"U"]), type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsdUt)

DFglsdUc= ur.ers(Canada[,"U"], type= "DF-GLS", model="const", lag.max=4)
summary(DFglsUc)

DFglsdUc= ur.ers(diff(Canada[,"U"]), type= "DF-GLS", model="const", lag.max=4)
summary(DFglsdUc)

###################################################################################
#Teste de quebra estrutural
raizquebra=ur.za((Canada[,"U"]), model='intercept', lag=3)
plot(raizquebra)
summary(raizquebra)
####################################################################################
###Vamos construir a equação do teste assim como um conjunto de variáveis sazonais
emprego = (Canada[,"e"])
emprego
length(emprego)

####################################################################################
# criando variável dummy trimestral
matr <- matrix(0, nr = 4, nc = 4)
matr[row(matr) == col(matr)] <- 1
matr
Djt <- rbind(matr,matr,matr,matr,matr,matr,matr,matr,matr,matr,matr,
             matr,matr,matr,matr,matr,matr,matr,matr,matr,matr); Djt
dim(Djt)


names(Djt)

dimnames(Djt) = list(c(),c("Q1", "Q2", "Q3", "Q4")) 
Djt
dim(Djt)

teste=data.frame(Djt)
attach(teste)
tri2= teste[1:83,2]
tri2
length(tri2)

tri3 = teste[1:83,3]
tri3

tri4 = teste[1:83,4]
Q4

demprego = (emprego[1:83])
demprego
length(demprego)

difemprego=diff(emprego)
difemprego
length(difemprego)

tempo = tempo <- seq(1:83) 
tempo

# teste de raíz unitária usando função dynlm com variáveis defasadas
modelo2 = dynlm(difemprego ~ demprego + tri2 +tri3 +tri4 + tempo + L(difemprego,1)+ L(difemprego,2)+ L(difemprego,3))
summary(modelo)

args(dynlm)

#################################################################################
#### Funções VAR ###########################################################
#causalidade de granger - refazer
require(MSBVAR)
require(lmtest)

names(Canada)
emprego=Canada$e
emprego

sareal = Canada$rw
sareal

desemprego = Canada$U
desemprego

produtividade = Canada$prod
produtividade 
#H0 não há causalidade de granger 
causalidade=grangertest(emprego ~ sareal, order=4)
causalidade
###################################################################################
### Obtendo os argumentos das funções 
args(VAR)
args(VARselect)
args(dynlm)
###################################################################################
## Determinação das defasagens do modelo a partir dos critérios de informação
VARselect(Canada, lag.max = 5, type = "const")

## Estimação do VAR (feita por MQO a partir de cada uma das equações) 
var_model = VAR(Canada, p=2, type = "const")
summary(var_model)

## se a constante não for significativa em nenhuma das equações, eliminamos esse termo

##análise gráfica
names(var_model)
plot(var_model) # apresenta um gráfico para cada variável

####################################################################################
##testes de diagnóstico

##normalidade dos resíduos

# verificando raízes do modelo estimado (fora do círculo unitário, valores menores do que 1)
roots(var_model)


#H0: residuo apresenta distribuição normal
normality.test(var_model, multivariate.only=TRUE)


## versão multivariada do teste de ljung-box (teste de portmanteau - estatistica Q)
args(serial)

##H0 não existe correlação até a ordem de defasagem n
serial.test(var_model, lags.pt=6, type="PT.adjusted")
serial.test(var_model, lags.pt=12, type="PT.adjusted")
serial.test(var_model, lags.pt=18, type="PT.adjusted")

##H0 teste individual não tem autocorrelação (teste LM)
serial.test(var_model, lags.bg=2, type="ES")
serial.test(var_model, lags.bg=4, tymobpe="ES")
serial.test(var_model, lags.bg=6, type="ES")


#Conditional heteroscedasticity
#H0 não há heterocedasticia
arch1_ca <- arch.test(var_model, lags.multi=3)
arch1_ca$arch.mul
arch1_ca <- arch.test(var_model, lags.multi=6)
arch1_ca$arch.mul
##################################################################################### 
##função de resposta a impulso
names(Canada)

impulso = irf(var_model, impulse = "rw", respomse = c("e","U"), boot = TRUE, cumulative= FALSE)
plot(impulso)

## - Decomposição da variância dos erros de predição
args(fevd)
decomp = fevd(var_model, n.ahead=10)
plot(decomp)
decomp  

# Modelo VAR com matriz triangular inferior
##################################################################################### 
## Para fazer previsão
prev = predict(var_model, n.ahead = 6, ci=0.95)
prev$fcst #coluna com o valor da previsão
plot(prev)

#######################################################################################
######################################################################################
# Cointegração Engle & Granger - duas séries - 17/07/23
#Dados Sorriso e Paranaguá
#Matriz Choleski
#   SO  PR
#SO 1  0
#PR 1  1

Sorriso = read.csv(file.choose())
Sorriso

attach(Sorriso)
summary(Sorriso)

is.numeric(SO)
is.numeric(PR)


lSO=log(SO)
lPR=log(PR)

length(lSO)
length(lPR)

#gerando dados do primeiro período
par(mfrow=c(1,1))

SO1=ts(lSO,start=1, end=800)
SO1
plot(SO1)

PR1=ts(lPR,start=1, end=800)
PR1
plot(PR1)

#gerando dados do segundo período
SO2=ts(lSO,start=801, end=3471)
SO2
plot(SO2)

PR2=ts(lPR,start=801, end=3471)
PR2
plot(PR2)

###testes de raiz unitária
#Série de preços Sorriso primeiro período 
adfSOt <- summary(ur.df(SO1, type="trend",  selectlags = c("AIC")))
adfSOt

adfSOc <- summary(ur.df(SO1, type="drift",  selectlags = c("AIC")))
adfSOc

adfSOn <- summary(ur.df(SO1, type="none",  selectlags = c("AIC")))
adfSOn

adfSOn <- summary(ur.df(diff(SO1), type="none",  selectlags = c("AIC")))
adfSOn

#Série de preços Paranaguá primeiro período 
adfPRt <- summary(ur.df(PR1, type="trend",  selectlags = c("AIC")))
adfPRt

adfPRc <- summary(ur.df(PR1, type="drift",  selectlags = c("AIC")))
adfPRc

adfPRn <- summary(ur.df(PR1, type="none",  selectlags = c("AIC")))
adfPRn

adfPRn <- summary(ur.df(diff(PR1), type="none",  selectlags = c("AIC")))
adfPRn


par(mfrow=c(1,2))
plot(SO1)
plot(PR1)

# ajustando o modelo de transmissão de preços - elasticidade de longo prazo
?lm
#Primeiro período
modelo1=lm(SO1 ~ PR1)
summary(modelo1)
anova(modelo1)

# modelo correção erro Sorriso PARANAGUÁ

modelo2 <- lm(SO2 ~ PR2)

ect2 = resid(modelo2)

plot(ect2)
adf.test(ect2)
length(ect2)

adfect2c <- summary(ur.df(ect2, type="drift",  selectlags = c("AIC")))
adfect2c

adfect2n <- summary(ur.df(ect2, type="none",  selectlags = c("AIC")))
adfect2n

#residuo vai rodar com as series na diff, perde uma observação

ectCE2<- resid(modelo1)[1:2670]
length(ectCE2)

dSO2 = diff(SO2)
dSO2
length(dSO2)

adf.test(dSO2)

dPR2 = diff(PR2)
dPR2
adf.test(dPR2)
length(dPR2)

adf.test(ect2)


dataCE <- cbind(dSO2, dPR2, ectCE2)


require (dynlm)
library(dynlm)

modeloCE2 = dynlm (dSO2 ~  dPR2 + L(ectCE2, 1) + L(dSO2, 1) + L(dPR2, 1) + L(dSO2, 2) + L(dPR2, 2) + L(dPR2, 3))
coef(modeloCE2)
summary(modeloCE2)

#######################################################################################
# Cointegração Johansen - 31/07/23
################# TESTE DE JOHANSEN ###############################################
#
#VECM e SVEC dados Canada
# 
data(Canada)
# Função ca.jo
# segue o mesmo exemplo, sendo a entrada "ecdet" aceita constante, tendencia ou nenhuma.

summary <- ca.jo(Canada, ecdet = "trend", type="trace", K=3, spec="transitory", season=4)
resultadoteste <- ca.jo(Canada, ecdet = "trend", type="trace", K=2, spec="transitory", season=4)
summary (resultadoteste)

vecm = ca.jo(Canada[, c("rw", "prod","e", "U")], type ="trace", ecdet ="trend",K=3, spec="transitory")
vecm.r1= cajorls(vecm, r=1)
summary(vecm.r1)


alphal=coef (vecm.r1$rlm) [1 , ] 
betal = vecm.r1$beta
resids=resid(vecm.r1$rlm)
N=nrow(resids)
sigmal=crossprod(resids)/N
vecm.r1@ZK
print(vecm.r1)
vecm.r1$beta
alphal


vecm = ca.jo(Canada[, c("rw", "prod","e", "U")], type ="trace", ecdet ="trend",K=3, spec="transitory")

LR<-matrix(NA,nrow=4,ncol=4)
LR[1,2:4]<-0
LR[2:4, 4]<-0
LR

SR=matrix(NA,nrow=4,ncol=4)
SR[4,2]<-0
SR

svec = SVEC(x=vecm,LR=LR,SR=SR,r=1,lrtest=FALSE,boot=TRUE, runs=100)
summary(svec)

data(finland)
# Função ca.jo
# segue o mesmos exemplos, sendo a entrada "ecdet" aceita constante, tendencia ou nenhuma.
sjf <- finland
sjf.vecm <- ca.jo(sjf, ecdet = "none", type="eigen", K=2, spec="longrun", season=4)
summary(sjf.vecm)

# Exemplo 2.
data(denmark)
sjd <- denmark[, c("LRM", "LRY", "IBO", "IDE")]
sjd.vecm <- ca.jo(sjd, ecdet = "const", type="trace", K=2, spec="longrun", season=4)
summary(sjd.vecm)

###############################################################################################################################################
#Teste de cointegração - Johansen dados Sorriso Paranaguá - aula 31/07/23
###############################################################################################################################################
#Período 1
y.matSO1PR1=data.frame(PR1, SO1)
vecmy.matSO1PR1= ca.jo(y.matSO1PR1, type="trace", ecdet="none", spec="longrun")
vecmy.matSO1PR1
summary(y.matSO1PR1)


#Período 2
y.matSO2PR2=data.frame(PR2, SO2)
vecmSO2PR2= ca.jo(y.matSO2PR2, type="trace", ecdet="none", spec="longrun")
vecmSO2PR2
summary(vecmSO2PR2)

#################################################################################################################################################
#VEC período 1
#################################################################################################################################################
dados1=cbind(PR1, SO1)
summary(dados1)
#Matriz Choleski
#   SO  PR
#SO  1  0
#PR  1  1


#Seleção defasagens modelo vec
VARselect(dados1, lag.max=15, type="both")

#teste de cointegração
y.matSO1PR1=data.frame(PR1, SO1)
vecmSO1PR1= ca.jo(y.matSO1PR1, type="trace", ecdet="none", spec="longrun",K=4)
summary(vecmSO1PR1)

#incorporando o vetor de correção de erro no VAR para estimar o VEC
veca=cajorls(vecmSO1PR1, r=1)
veca
summary(veca)
?cajorls

vec2var(vecmSO1PR1, r=1)
dev.off
plot(vecm1, names="lacr")


#VEC estrutural - SVEC
##Ajuste da matriz de relações contemporâneas 
LR<-matrix(NA,nrow=2,ncol=2)
LR
SR=NULL

vecmSO1PR1= ca.jo(y.matSO1PR1, type="trace", ecdet="const", spec="longrun",K=4)

SR=matrix(NA,nrow=2,ncol=2)
SR[2,1]<-0
SR

##Modelo
VECestrutural=SVEC(vecmSO1PR1,LR=LR,SR=SR,r=1,lrtest=FALSE,boot=FALSE)
summary(VECestrutural)
impulso=irf(VECestrutural,n.ahead=30,boot=TRUE)
impulso
decom=fevd(VECestrutural,n.ahead=30)
decom

#################################################################################################################################################
#VEC período 2
#################################################################################################################################################

dados2=cbind(SO2,PR2)
summary(dados2)

#Seleção defasagens modelo vec
VARselect(dados2, lag.max=15, type="both")

#teste de cointegração
y.matSO2PR2=data.frame(SO2, PR2)
vecmSO2PR2= ca.jo(y.matSO2PR2, type="trace", ecdet="const", spec="longrun",K=4)
summary(vecmSO2PR2)

#rodando o vec
veca2=cajorls(vecmSO2PR2, r=1)
veca2
summary(veca2)
?cajorls

vec2var2(vecmSO2PR2, r=1)
dev.off
plot(vecmSO2PR2, names="lacr")

#SVEC
##Ajuste da matriz de relações contemporâneas 
LR<-matrix(NA,nrow=2,ncol=2)
LR
SR=NULL
vecmSO2PR2= ca.jo(y.matSO2PR2, type="trace", ecdet="const", spec="longrun",K=4)
SR=matrix(NA,nrow=2,ncol=2)
SR[2,1]<-0
SR

##Modelo
VECestrutural2=SVEC(vecmSO2PR2,LR=LR,SR=SR,r=1,lrtest=FALSE,boot=FALSE)
summary(VECestrutural2)

impulso2=irf(VECestrutural2,n.ahead=30,boot=TRUE)
impulso2

decom2=fevd(VECestrutural2,n.ahead=30)
decom2
###########################################################################
#Dados Livea
#livea=read.csv(file.choose())
livea
attach(livea)
summary(livea)
str(livea) # verificar se estão como num (precisei colocar separador de unidade como. e decimal ,)

head(livea) #leitura das primeiras linhas da base.
names(livea) #lista das variáveis
dim(livea) #monstra a dimensão da base
livea[1:3, 1:4]  ## mostra linhas 1 a 4 e colunas 1 a 5 do data.frame

is.numeric(VPIB)
length(VPIB)
##############################################################################
#ANÁLISE GRÁFICA
##############################################################################
##Gráfico das séries em logaritmo
par(mfrow=c(2,2))
DPIB=ts(DPIB,frequency=4,start=c(1996,1))
plot.ts(DPIB)

DAGRO=ts(DAGRO,frequency=4,start=c(1996,1))
plot.ts(DAGRO)

AGTINT=ts(AGTINT,frequency=4,start=c(1996,1))
plot.ts(AGTINT)

CAMB=ts(CAMB,frequency=4,start=c(1996,1))
plot.ts(CAMB,frequency=4,start=c(1996,1))

JUR=ts(JUR,frequency=4,start=c(1996,1))
plot.ts(JUR,frequency=4,start=c(1996,1))

VAGRO=ts(VAGRO,frequency=4,start=c(1996,1))
plot.ts(VAGRO,frequency=4,start=c(1996,1))

VPIB=ts(VPIB,frequency=4,start=c(1996,1))
plot.ts(VPIB,frequency=4,start=c(1996,1))
##############################################################################
#transformação logaritmica das variáveis
##############################################################################
lDPIB=ts(log(DPIB),frequency=4,start=c(1996,1))
plot.ts(lDPIB)
acf(lDPIB)
pacf(lDPIB)
dlDPIB=diff(lDPIB)
acf(dlDPIB)

lDAGRO=ts(log(DAGRO),frequency=4,start=c(1996,1))
plot.ts(lDAGRO)
acf(lDAGRO)
pacf(lDAGRO)
dlDAGRO=diff(lDAGRO)
acf(dlDAGRO)
plot.ts(dlDAGRO)

lAGTINT=ts(log(AGTINT),frequency=4,start=c(1996,1))
plot.ts(lAGTINT)
acf(lAGTINT)
pacf(lAGTINT)
dlAGTINT=diff(lAGTINT)
acf(dlAGTINT)

lCAMB=ts(log(CAMB),frequency=4,start=c(1996,1))
plot.ts(lCAMB,frequency=4,start=c(1996,1))
acf(lCAMB)
pacf(lCAMB)
dlCAMB=diff(lCAMB)
acf(dlCAMB)

lJUR=ts(log(JUR),frequency=4,start=c(1996,1))
plot.ts(lJUR,frequency=4,start=c(1996,1))
acf(lJUR)
pacf(lJUR)
dlJUR=diff(lJUR)
acf(dlJUR)

lVAGRO=ts(log(VAGRO),frequency=4,start=c(1996,1))
plot.ts(lVAGRO,frequency=4,start=c(1996,1))
acf(lVAGRO)
pacf(lVAGRO)
dlVAGRO=diff(lVAGRO)
acf(dlVAGRO)

lVPIB=ts(log(VPIB),frequency=4,start=c(1996,1))
plot.ts(lVPIB,frequency=4,start=c(1996,1))
acf(lVPIB)
pacf(lVPIB)
dlVPIB=diff(lVPIB)
acf(dlVPIB)
##############################################################################
# teste de raíz unitária
################################################################################
#teste adf pacote urca

#SERIE DPIB
adfDPIBt <- summary(ur.df(lDPIB, type="trend", selectlags=c("AIC")))
adfDPIBt

adfDPIBc <- summary(ur.df(lDPIB, type="drift", selectlags=c("AIC")))
adfDPIBc

adfDPIBn <- summary(ur.df(lDPIB, type="none", selectlags=c("AIC")))
adfDPIBn

adfDPIBd <- summary(ur.df(dlDPIB, type="none", selectlags=c("AIC")))
adfDPIBd

#SERIE DAGRO
adfDAGROt <- summary(ur.df(lDAGRO, type="trend", selectlags=c("AIC")))
adfDAGROt

adfDAGROc <- summary(ur.df(lDAGRO, type="drift", lags= 4))
adfDAGROc

adfDAGROn <- summary(ur.df(lDAGRO, type="none", lags= 4))
adfDAGROn

adfDAGROd <- summary(ur.df(dlDAGRO, type="none", lags= 3))
adfDAGROd

?ur.df

#SERIE AGTINT
adfAGTINTt <- summary(ur.df(lAGTINT, type="trend", lags=1))
adfAGTINTt

adfAGTINTc <- summary(ur.df(lAGTINT, type="drift", lags=2))
adfAGTINTc

adfAGTINTn <- summary(ur.df(lAGTINT, type="none", lags=2))
adfAGTINTn

adfAGTINTd <- summary(ur.df(dlAGTINT, type="none", lags=1))
adfAGTINTd

#SERIE CAMB
adfCAMBt <- summary(ur.df(lCAMB, type="trend", lags=4))
adfCAMBt

adfCAMBc <- summary(ur.df(lCAMB, type="drift", lags=4))
adfCAMBc

adfCAMBn <- summary(ur.df(lCAMB, type="none", lags=1))
adfCAMBn

adfCAMBd <- summary(ur.df(dlCAMB, type="none", lags=1))
adfCAMBd

#SERIE JUR
adfJURt <- summary(ur.df(lJUR, type="trend", lags=3))
adfJURt

adfJURc <- summary(ur.df(lJUR, type="drift", lags=4))
adfJURc

adfJURn <- summary(ur.df(lJUR, type="none", lags=4))
adfJURn

adfJURd <- summary(ur.df(dlJUR, type="none", lags=3))
adfJURd

#SERIE VAGRO
adfVAGROt <- summary(ur.df(lVAGRO, type="trend", lags=4))
adfVAGROt

adfVAGROc <- summary(ur.df(lVAGRO, type="drift", lags=4))
adfVAGROc

adfVAGROn <- summary(ur.df(lVAGRO, type="none", lags=4))
adfVAGROn

adfVAGROd <- summary(ur.df(dlVAGRO, type="none", lags=3))
adfVAGROd

#SERIE VPIB
adfVPIBt <- summary(ur.df(lVPIB, type="trend", lags=4))
adfVPIBt

adfVPIBc <- summary(ur.df(lVPIB, type="drift", lags=4))
adfVPIBc

adfVPIBn <- summary(ur.df(lVPIB, type="none", lags=4))
adfVPIBn

adfVPIBd <- summary(ur.df(dlVPIB, type="none", lags=3))
adfVPIBd

# teste DF-GLS 
#SERIE DPIB
DFglDPIBt= ur.ers(lDPIB, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglDPIBt)

DFgldDPIBt= ur.ers(dlDPIB, type= "DF-GLS", model="trend", lag.max=4)
summary(DFgldDPIBt)

DFglDPIBc= ur.ers(lDPIB, type= "DF-GLS", model="const", lag.max=4)
summary(DFglDPIBc)

DFgdlDPIBc= ur.ers(dlDPIB, type= "DF-GLS", model="const", lag.max=4)
summary(DFgdlDPIBc)

#SERIE DAGRO
DFglsDAGROt= ur.ers(lDAGRO, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsDAGROt)

DFglsdDAGROt= ur.ers(dlDAGRO, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsdDAGROt)

DFglsDAGROc= ur.ers(lDAGRO, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsDAGROc)

DFglsdDAGROc= ur.ers(dlDAGRO, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsdDAGROc)

#SERIE AGTINT
DFglsAGTINTt= ur.ers(lAGTINT, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsAGTINTt)

DFglsdAGTINTt= ur.ers(dlAGTINT, type= "DF-GLS", model="trend", lag.max=3)
summary(DFglsdAGTINTt)

DFglsAGTINTc= ur.ers(lAGTINT, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsAGTINTc)

DFglsdAGTINTc= ur.ers(dlAGTINT, type= "DF-GLS", model="const", lag.max=3)
summary(DFglsdAGTINTc)

#SERIE CAMB
DFglsCAMBt= ur.ers(lCAMB, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsCAMBt)

DFglsdCAMBt= ur.ers(dlCAMB, type= "DF-GLS", model="trend", lag.max=3)
summary(DFglsdCAMBt)

DFglsCAMBc= ur.ers(lCAMB, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsCAMBc)

DFglsdCAMBc= ur.ers(dlCAMB, type= "DF-GLS", model="const", lag.max=3)
summary(DFglsdCAMBc)

#SERIE JUR
DFglsJURt= ur.ers(lJUR, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsJURt)

DFglsdJURt= ur.ers(dlJUR, type= "DF-GLS", model="trend", lag.max=3)
summary(DFglsdJURt)

DFglsJURc= ur.ers(lJUR, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsJURc)

DFglsdJURc= ur.ers(dlJUR, type= "DF-GLS", model="const", lag.max=3)
summary(DFglsdJURc)

#SERIE VAGRO
DFglsVAGROt= ur.ers(lVAGRO, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsVAGROt)

DFglsdVAGROt= ur.ers(dlVAGRO, type= "DF-GLS", model="trend", lag.max=3)
summary(DFglsdVAGROt)

DFglsVAGROc= ur.ers(lVAGRO, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsVAGROc)

DFglsdVAGROc= ur.ers(dlVAGRO, type= "DF-GLS", model="const", lag.max=3)
summary(DFglsdVAGROc)

#SERIE VPIB
DFglsVPIBt= ur.ers(lVPIB, type= "DF-GLS", model="trend", lag.max=4)
summary(DFglsVPIBt)

DFglsdVPIBt= ur.ers(dlVPIB, type= "DF-GLS", model="trend", lag.max=3)
summary(DFglsdVPIBt)

DFglsVPIBc= ur.ers(lVPIB, type= "DF-GLS", model="const", lag.max=4)
summary(DFglsVPIBc)

DFglsdVPIBc= ur.ers(dlVPIB, type= "DF-GLS", model="const", lag.max=3)
summary(DFglsdVPIBc)

#################################################################################################################################################
#VAR/VEC
#################################################################################################################################################
#################### Modelo VAR com as séries na diferença ####################
# análise de correlação cruzada DPIB e DAGRO
ccf(dlDPIB,dlDAGRO)
print(ccf(dlDPIB,dlDAGRO))

# análise de correlação cruzada DPIB e AGTINT
ccf(dlDPIB,dlAGTINT)
print(ccf(dlDPIB,dlAGTINT))

# análise de correlação cruzada DPIB e CAMB
ccf(dlDPIB,dlCAMB)
print(ccf(dlDPIB,dlCAMB))

###############################################################################
##  VAR
###############################################################################
# Construindo arquivo com dados no nível e nas diferemças
liveanivel=cbind(lCAMB,lAGTINT,lJUR,lDPIB,lDAGRO,lVAGRO,lVPIB)
liveanivel

dados = data.frame(lAGTINT,lCAMB,lJUR,lDPIB,lDAGRO,lVAGRO,lVPIB)
dados

dados_dif = data.frame(diff(dados$lAGTINT), diff(dados$lCAMB), diff(dados$lJUR),
                       diff(dados$lDPIB), diff(dados$lDAGRO), diff(dados$lVAGRO), diff(dados$lVPIB))
dados_dif

library(vars)
#Seleção defasagens modelo vec  - tentar rodar com vec2var
VARselect(dados_dif, lag.max=4, type="both")

#teste de cointegração
vecmlivea = ca.jo(dados, type="eigen", ecdet = "const", season = 4, spec = "longrun", K=3)
summary(vecmlivea)

vecm.r2=cajorls(vecmlivea , r = 2)
summary(vecm.r2)
vecm.r2

alphal=coef (vecm.r2$rlm) [1 , ] 
betal = vecm.r2$beta
resids=resid(vecm.r2$rlm)
N=nrow(resids)
sigmal=crossprod(resids)/N
vecm.r2@ZK
print(vecm.r2)
vecm.r2$beta

## estatística t para alpha
#obs: começar a rodar o código do alga do parenteses central pro resto - tem que definir vecm@ZK%*%beta primeiro!
alphal.se=sqrt(solve(crossprod(cbind(vecmlivea@ZK%*%betal,vecmlivea@Z1) ) )[1,1]*diag(sigmal)) 
alphal.t=alphal/alphal.se
alphal.t

## estatística t para beta

betal.se= sqrt(diag(kronecker(solve(crossprod(vecmlivea@RK[, -1])),solve(t(alphal)%*%solve(sigmal)%*% alphal))))
betal1=betal [,1] #coluna 1 de beta
betal2=betal[,2] #coluna 2 de beta
betal1
betal2
betal1[-1]#retira o coeficiente normalizado do primeiro vetor
betal2[-2]#retira o coeficiente normalizado do segundo vetor
betal.t1 = c(NA,betal1[-1]/betal.se)
betal.t1
betal.t2 = c(NA,betal2[-2]/betal.se)
betal.t2


#testes de diagnóstico
mod_diag = vec2var (vecmlivea, r=2)
mod_diag


#h0: residuo apresenta distribuição normal
normality.test(mod_diag, multivariate.only=TRUE)

#testes de autocorrelação
#H0 não existe correlação até a ordem de defasagem n
serial.test(mod_diag, lags.pt = 6, type="PT.adjusted")
serial.test(mod_diag, lags.pt = 8, type="PT.adjusted")
serial.test(mod_diag, lags.pt = 12, type="PT.adjusted")

serial.test(mod_diag, lags.bg = 1, type="ES")
serial.test(mod_diag, lags.bg = 2, type="ES")
serial.test(mod_diag, lags.bg = 3, type="ES")

#Conditional heteroscedasticity
#H0 não há heterocedasticia
arch1_br <- arch.test(mod_diag, lags.multi=3)
arch1_br$arch.mul
arch1_br <- arch.test(mod_diag, lags.multi=6)
arch1_br$arch.mul

############################################################################################
# Estimando o modelo incluindo o vetor de cointegração como variável exógena  e estimando a 
# elasticidade de curto prazo a partir da matriz de relações contemporâneas

# Construindo arquivo com dados no nível e nas diferemças
liveanivel=cbind(lCAMB,lAGTINT,lJUR,lDPIB,lDAGRO,lVAGRO,lVPIB)
liveanivel

dados = data.frame(lCAMB,lAGTINT,lJUR,lDPIB,lDAGRO,lVAGRO,lVPIB)
dados

dados_dif = data.frame(diff(dados$lCAMB), diff(dados$lAGTINT), diff(dados$lJUR),
                       diff(dados$lDPIB), diff(dados$lDAGRO), diff(dados$lVAGRO), diff(dados$lVPIB))
dados_dif

#identificando o número de defasagens 
VARselect(dados_dif, lag.max = 4, type="const", season = 4)

#Testes de cointegração
#Estimando o vetor de cointegração
teste_max = ca.jo(dados, type="eigen", ecdet = "const", K=3,  season = 4, spec = "transitory")
summary(teste_max)

#ajustando o modelo VEC
mod_vec = cajorls(teste_max , r=2)
print(mod_vec)
mod_vec$beta

# Construindo os vetores de cointegração ect1 e ect 2 
vetor_coint1 = (lag(dados$lCAMB))*mod_vec$beta[1,1] + mod_vec$beta[2,1]*lag(dados$lAGTINT)+mod_vec$beta[3,1]*lag(dados$lJUR)+
  mod_vec$beta[4,1]*lag(dados$lDPIB)+mod_vec$beta[5,1]*lag(dados$lDAGRO)+mod_vec$beta[6,1]*lag(dados$lVAGRO)+
  mod_vec$beta[7,1]*lag(dados$lVPIB)+ mod_vec$beta[8,1]
vetor_coint1
mod_vec$beta[7,1]

vetor_coint2 = (lag(dados$lCAMB))*mod_vec$beta[1,2] + mod_vec$beta[2,2]*lag(dados$lAGTINT)+mod_vec$beta[3,2]*lag(dados$lJUR)+
  mod_vec$beta[4,2]*lag(dados$lDPIB)+mod_vec$beta[5,2]*lag(dados$lDAGRO)+mod_vec$beta[6,2]*lag(dados$lVAGRO)+
  mod_vec$beta[7,2]*lag(dados$lVPIB)+mod_vec$beta[8,2]
vetor_coint2


length(vetor_coint2)

cointegra = data.frame(vetor_coint1[1:89],vetor_coint2[1:89])
cointegra

var_coint = VAR (dados_dif, p=3, type="none", season = 4, exogen=cointegra[1:89,])
var_coint
print(var_coint)

matriz.a = diag(7)
matriz.a[2,1]=NA
matriz.a[3,1]=NA
matriz.a[3,2]=NA
matriz.a[4,1]=NA
matriz.a[4,2]=NA
matriz.a[4,3]=NA
matriz.a[5,1]=NA
matriz.a[5,2]=NA
matriz.a[5,3]=NA
matriz.a[5,4]=NA
matriz.a[6,1]=NA
matriz.a[6,2]=NA
matriz.a[6,3]=NA
matriz.a[6,4]=NA
matriz.a[6,5]=NA
matriz.a[7,1]=NA
matriz.a[7,2]=NA
matriz.a[7,3]=NA
matriz.a[7,4]=NA
matriz.a[7,5]=NA
matriz.a[7,6]=NA
matriz.a

args(optim)

svar_coint = SVAR (var_coint, Amat=matriz.a, Bmat=NULL, estmethod = "scoring", method = "BFGS", hessian=TRUE, max.iter = 100, conv.crit = 0.1e-10)
svar_coint$A
summary(svar_coint)
print(svar_coint)

impulso=irf(svar_coint,n.ahead = 4,boot=TRUE)
impulso
plot(impulso)

decom=fevd(VECestrutural,n.ahead = 4)
decom
plot(decom)