A seguir, encotram-se todos os códigos utilizados para a confecção do artigo " A dinãmica do poder de compra através do salário mínimo real: uma análise temporal de 2000 a 2019"
# PACOTES ------------------------------------------------------------------------------------------------------------------------------------------------------
#install.packages("magic")
library(readxl)
library(magic)
# DADOS --------------------------------------------------------------------------------------------------------------------------------------------------------
setwd("E:\\dinamico")
sam <- read_excel("Salário Mínimo Real.xls")
y = sam[715:nrow(sam),2]
y <- (y$`Salário mínimo real - R$ - Instituto de Pesquisa Econômica Aplicada(IPEA) - GAC12_SALMINRE12 -`)
# Trabalhando de 2000.1 a 2019.10
# VISUALIZAÇÃO -------------------------------------------------------------------------------------------------------------------------------------------------
#traço da série
ts.plot(y,col=2)
#autocorrelação
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,2,1)+0.1)
acf(y,lag.max = 120,ylab="Autocorrelação",bty="n",main="",xaxt='n',xlab="",cex.lab = 3,cex.axis=3, cex.main=3, cex.sub=3,lwd=3) #não estacionário
axis(1, at=seq(0,120,by=10), las=2,cex.axis=3)
title(xlab = "Lag", cex.lab = 3,line = 6)
adif <- diff(as.vector(y),lag = 1) #fazendo diferenciação para remover estacionariedade
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,2,1)+0.1)
acf(adif,lag.max = 120,ylab="Autocorrelação",bty="n",main="",xaxt='n',xlab="",cex.lab = 3,cex.axis=3, cex.main=3, cex.sub=3,lwd=3) #não estacionário
axis(1, at=seq(0,120,by=10), las=2,cex.axis=3)
title(xlab = "Lag", cex.lab = 3,line = 6)
#harmônico
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,2,1)+0.1)
plot(y,col=2,main="",xaxt='n',yaxt='n',xlab="",ylab="",cex.lab = 3,cex.axis=3, cex.main=3, cex.sub=3,lwd=3,type="l",bty="n")
abline(v=seq(1,length(y)+1,by = 12),lty=2,lwd=2)
axis(1, at=seq(0,240,by=12), las=2,cex.axis=3)
axis(2, at=seq(500,1000,by=100), las=1,cex.axis=3)
title(xlab = "Tempo",ylab = "y", cex.lab = 3,line = 7)
#o comportamento não parece ser harmônico
#porém o mais próximo desse comportamento é o 1º
#removendo o ultimos 13 meses para previsão
h <- 0
l <- length(y)
yt <- as.vector(y[1:(l-h)])
y2 <- as.vector(y[(l-h+1):l])
p <- 12
G1 <- matrix(c(1,0,1,1),2)
w <- 2*pi/p
G2 <- matrix(c(cos(w),-sin(w),sin(w),cos(w)),2)
G <- adiag(G1,G2)
# FILTRAGEM ----
j <- 4
N <- length(yt)
m0 <- matrix(rep(0,j),ncol=1)
C0 <- diag(100,j)
delta <- 0.95
n <- 1
d0 <- 1
K0 <- d0/n
Ft <- t(cbind(rep(1,l),0,1,0))
#t(cbind(rep(1,N),0,1,0))
at <- matrix(NA,ncol=N,nrow=j)
mt <- matrix(NA,ncol=N,nrow=j)
A <- matrix(NA,ncol=N,nrow=j)
qt <- rep(NA,N)
ft <- rep(NA,N)
et <- rep(NA,N)
d <- rep(NA,N)
rt <- list()
Ct <- list()
Kt <- rep(NA,N)
for(i in 1:N){
if(i==1){
at[,i] <- G%*%m0
rt[[i]] <- (G%*%C0%*%t(G))/delta
ft[i] <- t(Ft[,i]) %*% as.matrix(at[,i])
qt[i] <- (t(Ft[,i])%*%rt[[i]]%*%Ft[,i])+K0
A[,i] <- rt[[i]]%*%Ft[,i]*as.numeric(1/qt[i])
et[i] <- yt[i]-ft[i]
mt[,i] <- at[,i] +as.matrix(A[,i])%*%et[i]
d[i] <- d0+K0*(et[i]^2)*(1/qt[i])
n <- n+1
Kt[i] <- d[i]/n
Ct[[i]] <- (Kt[i]/K0)*(rt[[i]] - as.matrix(A[,i])%*%t(as.matrix(A[,i]))*as.numeric(qt[i]))
}else{
at[,i] <- G%*%mt[,i-1]
rt[[i]] <- (G%*%Ct[[(i-1)]]%*%t(G))/delta
ft[i] <- t(Ft[,i]) %*% as.matrix(at[,i])
qt[i] <- (t(Ft[,i])%*%rt[[i]]%*%Ft[,i])+Kt[(i-1)]
A[,i] <- rt[[i]]%*%Ft[,i]*as.numeric(1/qt[i])
et[i] <- yt[i]-ft[i]
mt[,i] <- at[,i] +as.matrix(A[,i])%*%et[i]
d[i] <- d[(i-1)]+Kt[(i-1)]*(et[i]^2)*(1/qt[i])
n <- n+1
Kt[i] <- d[i]/n
Ct[[i]] <- (Kt[i]/Kt[(i-1)])*(rt[[i]] - as.matrix(A[,i])%*%t(as.matrix(A[,i]))*as.numeric(qt[i]))
}
}
cred <- cbind(ft+qnorm(0.025)*sqrt(qt),ft+qnorm(0.975)*sqrt(qt))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,8,2)+0.1)
plot(ft, bty="l",type="l", lwd=3, axes = FALSE, ann = FALSE,lty=2,ylim=c(0,1100))
lines(yt, col=rgb(0.1,0.8,0.1,1) , lwd=3 , pch=18 , type="l")
lines(cred[,1], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
lines(cred[,2], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
title(ylab = "y", cex.lab = 2.5,line = 7)
title(xlab = "Tempo", cex.lab = 2.5,line = 7)
axis(1, at=seq(0,length(y),by=7), las=2,cex.axis=3)
axis(2, at=seq(0,1000,by=100), las=1,cex.axis=3)
legend("bottomleft",
legend = c("Previsão on-line", "Valor Verdadeiro","95% IC"),
col = c("black",
rgb(0.1,0.8,0.1,1) ,
rgb(0.8,0.2,0.2,1) ,
rgb(0.8,0.2,0.2,1)),
bty = "n",
lty = c(2,1,2),
pt.cex = 2,
cex = 2.5,
lwd=3,
text.col = "black",
horiz = F ,
inset = c(0.6, 0.2))
# SUAVIZAÇÃO ----
bt <- list()
st <- matrix(NA,nrow=j,ncol=N)
st[,N] <- mt[,N]
St <- list()
St[[N]] <- Ct[[N]]
for(i in (N-1):1){
svd_R <- svd(rt[[i+1]])
R_inv <- svd_R$v%*%diag(1/svd_R$d)%*%t(svd_R$u)
bt[[i]] <- Ct[[i]]%*%t(G)%*%R_inv
st[,i] <- mt[,i] + bt[[i]]%*%(st[,i+1]-at[,(i+1)])
St[[i]] <- Ct[[i]] - bt[[i]] %*% (rt[[i+1]] - St[[i+1]]) %*% t(bt[[i]])
}
S1 <- sapply(St, function(x){x[1,1]})
cred <- cbind(st[1,]+qnorm(0.025)*sqrt(S1),st[1,]+qnorm(0.975)*sqrt(S1))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,8,2)+0.1)
plot(st[1,], type="l", lwd=3, axes = FALSE, ann = FALSE,ylim=c(-100,1100))
lines(cred[,1], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
lines(cred[,2], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
abline(h=0,lty=2,lwd=3)
title(ylab = expression(mu), cex.lab = 2.5,line = 7)
title(xlab = "Tempo", cex.lab = 2.5,line = 7)
axis(1, at=seq(0,length(y),by=12), las=2,cex.axis=3)
axis(2, at=seq(0,1000,by=100), las=1,cex.axis=3)
legend("bottomleft",
legend = c("Estimativa suavizada","95% IC","Zero"),
col = c("black",rgb(0.8,0.2,0.2,1),"black"),
bty = "n",
lty = c(1,2,2),
pt.cex = 1,
cex = 2,
lwd=3,
seg.len=1,
text.col = "black",
horiz = F ,
inset = c(0.55, 0.2))
S2 <- sapply(St, function(x){x[2,2]})
cred <- cbind(st[2,]+qnorm(0.025)*sqrt(S2),st[2,]+qnorm(0.975)*sqrt(S2))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,8,2)+0.1)
plot(st[2,], type="l", lwd=3, axes = FALSE, ann = FALSE,ylim=c(-3,30))
lines(cred[,1], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
lines(cred[,2], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
abline(h=0,lty=2,lwd=3)
title(ylab = expression(beta), cex.lab = 2.5,line = 7)
title(xlab = "Tempo", cex.lab = 2.5,line = 7)
axis(1, at=seq(0,length(y),by=12), las=2,cex.axis=3)
axis(2, at=seq(-3,30,by=3), las=1,cex.axis=3)
legend("bottomleft",
legend = c("Estimativa suavizada","95% IC","Zero"),
col = c("black",rgb(0.8,0.2,0.2,1),"black"),
bty = "n",
lty = c(1,2,2),
pt.cex = 2,
cex = 2,
lwd=3,
seg.len=1,
text.col = "black",
horiz = F ,
inset = c(0.55, 0.2))
S3 <- sapply(St, function(x){x[3,3]})
cred <- cbind(st[3,]+qnorm(0.025)*sqrt(S3),st[3,]+qnorm(0.975)*sqrt(S3))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,8,2)+0.1)
plot(st[3,], type="l", lwd=3, axes = FALSE, ann = FALSE,ylim=c(-80,40))
lines(cred[,1], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
lines(cred[,2], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
abline(h=0,lty=2,lwd=3)
title(ylab = expression(gamma), cex.lab = 2.5,line = 7)
title(xlab = "Tempo", cex.lab = 2.5,line = 7)
axis(1, at=seq(0,length(y),by=12), las=2,cex.axis=3)
axis(2, at=seq(-80,40,by=10), las=1,cex.axis=3)
legend("bottomleft",
legend = c("Estimativa suavizada","95% IC","Zero"),
col = c("black",rgb(0.8,0.2,0.2,1),"black"),
bty = "n",
lty = c(1,2,2),
pt.cex = 2,
cex = 2,
lwd=3,
seg.len=1,
text.col = "black",
horiz = F ,
inset = c(0.55, 0.05))
S4 <- sapply(St, function(x){x[4,4]})
cred <- cbind(st[4,]+qnorm(0.025)*sqrt(S4),st[4,]+qnorm(0.975)*sqrt(S4))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,9,8,2)+0.1)
plot(st[4,], type="l", lwd=3, axes = FALSE, ann = FALSE,ylim=c(-80,40))
lines(cred[,1], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
lines(cred[,2], col=rgb(0.8,0.1,0.1,1) , lwd=3 , pch=18 , type="b")
abline(h=0,lty=2,lwd=3)
title(ylab = expression(delta), cex.lab = 2.5,line = 7)
title(xlab = "Tempo", cex.lab = 2.5,line = 7)
axis(1, at=seq(0,length(y),by=12), las=2,cex.axis=3)
axis(2, at=seq(-80,40,by=10), las=1,cex.axis=3)
legend("bottomleft",
legend = c("Estimativa suavizada","95% IC","Zero"),
col = c("black",rgb(0.8,0.2,0.2,1),"black"),
bty = "n",
lty = c(1,2,2),
pt.cex = 2,
cex = 2,
lwd=3,
seg.len=1,
text.col = "black",
horiz = F ,
inset = c(0.55, 0.05))
#previsão
h = 12
N = length(y) - h
y2 <- as.vector(y[(l-h+1):l])
a_prev = matrix(NA, h, nrow(Ft))
f_prev = vector()
R_prev = list()
Q_prev = vector()
theta_prev = matrix(NA, h, nrow(Ft))
a_prev[1,] = mt[,N]
R_prev[[1]] = Ct[[N]]
f_prev[1] <- t(Ft[,N+1]) %*% a_prev[1, ]
Q_prev[1] <- t(Ft[,N+1]) %*% R_prev[[1]] %*% Ft[,N+1] + Kt[[N]]
for (t in 2:h){
a_prev[t, ] <- G %*% a_prev[t-1, ]
R_prev[[t]] <- G %*% R_prev[[t-1]] %*% t(G)/delta
f_prev[t] <- t(Ft[,N+t]) %*% a_prev[t, ]
Q_prev[t] <- t(Ft[,N+t]) %*% R_prev[[t]] %*% Ft[,N+t] + Kt[[N]]
}
cred <- cbind(f_prev+qnorm(0.025)*sqrt(Q_prev),f_prev+qnorm(0.975)*sqrt(Q_prev))
# PREVISAO ----
x11()
par(mgp=c(5,1,0))
par(mar=c(9,10,2,1)+0.1)
plot(yt,col=rgb(0.1,0.8,0.1,1),main="",xaxt='n',yaxt='n',xlab="",ylab="",cex.lab = 3,cex.axis=3, cex.main=3, cex.sub=3,lwd=3,type="l",bty="n",xlim = c(0,250),ylim=c(500,1100))
abline(v=N,lty=2,lwd=3)
lines((N+1):(N+h),f_prev, col=rgb(0.1,0.1,0.8,1),lwd=3,type="b",bty="l",pch=20,cex=1)
lines((N+1):(N+h),cred[,1], col=rgb(0.8,0.1,0.1,1),lwd=3,type="b",bty="l",pch=18,cex=1,lty=3)
lines((N+1):(N+h),cred[,2], col=rgb(0.8,0.1,0.1,1),lwd=3,type="b",bty="l",pch=18,cex=1,lty=3)
axis(1, at=seq(0,250,by=10), las=2,cex.axis=3)
axis(2, at=seq(500,1100,by=100), las=1,cex.axis=3)
title(xlab = "Tempo",ylab = "y", cex.lab = 3,line = 7)
legend("bottomleft",
legend = c("Corte da previsão","Valor previsto", "Valor Verdadeiro","95% IC"),
col = c("black",
rgb(0.1,0.8,0.1,1) ,
rgb(0.1,0.1,0.8,1),
rgb(0.8,0.2,0.2,1) ,
rgb(0.8,0.2,0.2,1)),
bty = "n",
lty = c(1,1,2,2),
pch = c(NA,NA,20,18),
pt.cex = 2,
cex = 2.5,
lwd=3,
text.col = "black",
horiz = F ,
inset = c(0.01, 0.5))
x11()
par(mgp=c(5,1,0))
par(mar=c(9,10,2,1)+0.1)
plot(yt[N:(N+h-1)],type="b",bty="l",pch=4,cex=3,col=rgb(0.1,0.8,0.1,1),main="",xaxt='n',yaxt='n',xlab="",ylab="",lwd=3,bty="n",ylim=c(min(cred),max(cred)),xlim=c(1,20))
abline(v=N,lty=2,lwd=3)
lines(f_prev, col=rgb(0.1,0.1,0.8,1), lwd=3,type="b",bty="l",pch=20,cex=3)
lines(cred[,1], col=rgb(0.8,0.1,0.1,1), lwd=3,type="b",bty="l",pch=18,cex=3,lty=3)
lines(cred[,2], col=rgb(0.8,0.1,0.1,1), lwd=3,type="b",bty="l",pch=18,cex=3,lty=3)
axis(1, at=seq(1,12,by=1),labels=seq(1,12,by=1), las=2,cex.axis=3)
axis(2, at=seq(940,1080,by=10), las=1,cex.axis=3)
title(xlab = "Tempo",cex.lab = 3,line = 7)
title(ylab = "y", cex.lab = 3,line = 8)
legend("bottomleft",
legend = c("Valor previsto", "Valor Verdadeiro","95% IC"),
col = c(rgb(0.1,0.8,0.1,1) ,
rgb(0.1,0.1,0.8,1),
rgb(0.8,0.2,0.2,1) ,
rgb(0.8,0.2,0.2,1)),
bty = "n",
lty = c(1,2,2),
pch = c(4,20,18,18),
seg.len=0.7,
cex = 2.5,
lwd=3,
text.col = "black",
horiz = F ,
inset = c(0.6, 0.3))
#Erro Percentual Absoluto Médio (MAPE)
(MAPE <- mean(abs(y2-f_prev)/abs(y2))*100)
#Erro Percentual Absoluto Médio (MAPE)
(MAPE2 <- mean(abs(y-ft)/abs(y))*100)
#salvar com um nomeqlqr.tex(da pra abrir com o bloco de notas pra pegar a table ja pronta)
print(xtable(a, type = "latex"), file = file.choose())