Modelos Dinâmicos Lineares (Graduação) - Lab 1
Introdução
O exercício a seguir é sobre o modelo polinomial de primeira ordem.
- Gere séries temporais para diferentes valores de \(V_t = V\) e \(W_t = W\), \(\forall t\). Faça o gráfico da série do nível \(\mu_t\) e das observações \(Y_t\) para cada caso considerado acima e compare o comportamento dos gráficos considerando a relação entre \(V\) e \(W\).
- Assumindo \(V=1,1\) e \(W=0,5\) gere uma série de tamanho T = 123 e em seguida ajuste o modelo utilizando o procedimento de inferência squencial
- Baseado no modelo polinomial de primeira ordem constante, verifique suas propriedades assintóticas.
1) Relação sinal-ruído
Tabela 1 - Valores para V e W investigados
| V | W | Razão |
|---|---|---|
| 1 | 9 | 9.0000000 |
| 2 | 8 | 4.0000000 |
| 3 | 7 | 2.3333333 |
| 4 | 6 | 1.5000000 |
| 5 | 5 | 1.0000000 |
| 6 | 4 | 0.6666667 |
| 7 | 3 | 0.4285714 |
| 8 | 2 | 0.2500000 |
| 9 | 1 | 0.1111111 |
Para avaliar o comportamento do modelo polinomial de primeira ordem, propõe-se um estudo simulado a partir da variação da quantidade \(R = \cfrac{W}{V}\). Para fins de simplificação de notação, este R será chamado de Razão, a doravante. Inicialmente, é de interesse avaliar as mudanças no comportamento de \(Y_t\) para quando a razão está em um raio crescendo de 0.1 até 9 e obter evidências visuais do comportamento da relação sinal-ruído.
Gráfico 1 - Evolução da Razão
O Gráfico 1 apresenta a relação visual entre as observações e o seu respectivo nível. Verificou-se que para razões altas, isto é, um \(W\) maior que o \(V\), as observações são possivelmente “dominadas” pela tendência. Analogamente, para uma relação sinal-ruído baixa, a viariabilidade de \(Y_t\) parece decorrer principalmente dos erros aleatórios com variância \(V_t\).
2) Ajuste e previsão
Gráfico 2 - Série observada vs Série ajustada
O Gráfico 2 representa a comparação visual entre a série observada e a série ajustada, para 120 observações. As linhas contínuas em negrito representam o Intervalo de Credibilidade para a Distribuição Preditiva, considerando 95% de probabiliadade. Verificou-se que a série ajustada é adequada ao acompanhar o comportamento da série verdadeira e manter-se dentro dos envelopes formados pelo intervalo de credibilidade.
Gráfico 3 - Previsão da Série ajustada com horizonte de 3 passos
o Gráfico 3 é um recorte do Gráfico 2, para as observações compreendidas entre o tempo 120 e 123. Tem-se por objetivo aplicar um zoom nas observações para que esta tenha uma visualização mais simples. Observou-se que a previsão para um horizonte de tempo igual 3 é constante e igual a \(m_{120}\), como definido no modelo.
3) Visualização assintótica
Gráfico 4 - \(A_t \xrightarrow{} A\)
Gráfico 5 - \(C_t \xrightarrow{} C\)
Gráfico 6 - \(R_t \xrightarrow{} R\)
Gráfico 7 - \(Q_t \xrightarrow{} Q\)
Os Gráficos 4, 5, 6 e 7 representam avaliações empíricas a respeito das propriedades de convergência assintótica do modelo constante. Verificou que para todos os parâmetros investigados, o comportamento esperado é atendido. Ademais, a convergência é razoavelmente próxima (em termos visuais) para valores de \(t>10\), em todos os casos observados.
Gráfico 8 - Previsão da Série ajustada com horizonte de 3 passos
Por fim, Gráfico 8 é uma representação visual das previsões para o tempo 121, 122 e 123. Verificou-se que as observações verdadeiras encontram-se dentro do Intervalo de Credibilidade produzido a partir da Distribuição Preditiva, com 95% de probabilidade, sendo este um indicativo de adequabilidade.
Apêndice A - Códigos
1)
set.seed(07082019)
MDPO = function(n,V,W){
y = vector()
mu = vector()
mu[1] = rnorm(1,0,10)
y[1] = mu[1] + rnorm(1,0,W)
for (t in 2:n){
mu[t] = mu[t-1] + rnorm(1,0,W)
y[t] = mu[t] + rnorm(1,0,V)
}
return(cbind(y,mu))
}
v = 1:9
w = 9:1
vw = cbind(v,w,w/v)
colnames(vw) = c("V","W","Razão")
kable(vw)%>%
kable_styling(bootstrap_options = c("striped", "hover"))
p1 = MDPO(100,v[1],w[1])
p2 = MDPO(100,v[2],w[2])
p3 = MDPO(100,v[3],w[3])
p4 = MDPO(100,v[4],w[4])
p5 = MDPO(100,v[5],w[5])
p6 = MDPO(100,v[6],w[6])
p7 = MDPO(100,v[7],w[7])
p8 = MDPO(100,v[8],w[8])
p9 = MDPO(100,v[9],w[9])
par(mfrow=c(3,1))
plot(p1[,1],type="l",main=round(vw[1,3],2),xlab="",ylab="")
lines(p1[,2],type="l",main=round(vw[1,3],2),col="red",xlab="",ylab="")
plot(p2[,1],type="l",main=round(vw[2,3],2),xlab="",ylab="")
lines(p2[,2],type="l",main=round(vw[2,3],2),col="red",xlab="",ylab="")
plot(p3[,1],type="l",main=round(vw[3,3],2),xlab="",ylab="")
lines(p3[,2],type="l",main=round(vw[3,3],2),col="red",xlab="",ylab="")
par(mfrow=c(3,1))
plot(p4[,1],type="l",main=round(vw[4,3],2),xlab="",ylab="")
lines(p4[,2],type="l",main=round(vw[4,3],2),col="red",xlab="",ylab="")
plot(p5[,1],type="l",main=round(vw[5,3],2),xlab="",ylab="")
lines(p5[,2],type="l",main=round(vw[5,3],2),col="red",xlab="",ylab="")
plot(p6[,1],type="l",main=round(vw[6,3],2),xlab="",ylab="")
lines(p6[,2],type="l",main=round(vw[6,3],2),col="red",xlab="",ylab="")
par(mfrow=c(3,1))
plot(p7[,1],type="l",main=round(vw[7,3],2),xlab="",ylab="")
lines(p7[,2],type="l",main=round(vw[7,3],2),col="red",xlab="",ylab="")
plot(p8[,1],type="l",main=round(vw[8,3],2),xlab="",ylab="")
lines(p8[,2],type="l",main=round(vw[8,3],2),col="red",xlab="",ylab="")
plot(p9[,1],type="l",main=round(vw[9,3],2),xlab="",ylab="")
lines(p9[,2],type="l",main=round(vw[9,3],2),col="red",xlab="",ylab="")2)
set.seed(42)
s = MDPO(123,1.1,0.5)
y = s[1:120,1]
V = 1.1
W = 0.5
m0 = 0
C0 = 100
m = vector()
R = vector()
Q = vector()
C = vector()
f = vector()
A = vector()
e = vector()
cred = matrix(NA,120,2)
f[1] = m0
R[1] = C0 + W
Q[1] = R[1] + V
A[1] = R[1]/Q[1]
e[1] = y[1] - f[1]
C[1] = A[1]*V
m[1] = m0 + A[1]*e[1]
#cred[1,] = qnorm(c(0.025,0.975),f[1],sqrt(Q[1]))
for (t in 2:120){
f[t] = m[t-1]
R[t] = C[t-1] + W
Q[t] = R[t] + V
A[t] = R[t]/Q[t]
e[t] = y[t] - f[t]
C[t] = A[t]*V
m[t] = m[t-1] + A[t]*e[t]
cred[t,] = qnorm(c(0.025,0.975),f[t],sqrt(Q[t]))
}
res = s[121:123,1] - m[120]
prev = m[120]
Ct = C[120]
Vt = V
ic1 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*1))
ic2 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*2))
ic3 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*3))
ic = rbind(ic1,ic2,ic3)
cred = rbind(cred,ic)
m[121:123] = prev
plot(s[,1],type="l",main="",xlab="x",col="blue",ylim=c(3,20),xlim=c(0,125),ylab="")
lines(m,type="l",col="red",lwd=2,lty=2)
lines(cred[,1],type="l",lwd=2,lty=1)
lines(cred[,2],type="l",lwd=2,lty=1)
abline(v=120)
legend(90, 18, legend=c("Observado", "Ajustado","IC"),col=c("blue", "red","black"), lty=c(1,2,1),cex=0.8)
plot(s[,1],type="l",main="",xlab="x",col="blue",ylim=c(3,15),xlim=c(100,125),ylab="")
lines(m,type="l",col="red",lwd=2,lty=2)
lines(cred[,1],type="l",lwd=2,lty=1)
lines(cred[,2],type="l",lwd=2,lty=1)
abline(v=120)
legend(120.5, 15, legend=c("Observado", "Ajustado","IC"),col=c("blue", "red","black"), lty=c(1,2,1),cex=0.8)3)
r = W/V
A_conv = r/2*(sqrt(1+4/r)-1)
plot(A,type="l",ylim=c(0,4))
abline(A_conv,0,col=2,lty=2)
legend(100, 4, legend=c(expression(At),expression(A )),col=c( "black","red"), lty=c(1,2),cex=0.8)plot(C,type="l",ylim=c(0,4))
abline(A_conv*V,0,col=2,lty=2)
legend(100, 4, legend=c(expression(Ct),expression(C )),col=c( "black","red"), lty=c(1,2),cex=0.8)plot(R,type="l",ylim=c(0,4))
abline(A_conv*V + W,0,col=2,lty=2)
R_conv = A_conv*V + W
legend(100, 4, legend=c(expression(Rt),expression(R)),col=c( "black","red"), lty=c(1,2),cex=0.8)plot(Q,type="l",ylim=c(0,4))
abline(R_conv + V,0,col=2,lty=2)
legend(100, 4, legend=c(expression(Qt),expression(Q)),col=c( "black","red"), lty=c(1,2),cex=0.8)res = s[121:123,1] - m[120]
prev = m[120]
Ct = C[120]
Vt = V
ic1 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*1))
ic2 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*2))
ic3 = qnorm(c(0.025,0.975),prev,sqrt(Ct+Vt+W*3))
ic = rbind(ic1,ic2,ic3)
plot(ic[,1],ylim=c(0.9,14),xlim=c(1,4),pch=25,ylab=expression(Yt))
points(ic[,2],pch=24)
points(1:3,rep(prev,3),pch="x",col="blue")
points(1:3,s[121:123,1],pch=18,col=2)
legend(3.4, 10, legend=c("Upper", "Previsto","Observado","Lower"),col=c("black", "blue","red","black"), pch=c(24,4,18,25),cex=0.8)