###########################################################################
#
# ECONOMETRIA AVANCADA 
# INSPER
# QUARTA LISTA DE EXERCICIOS
# PRIMEIRO SEMESTRE 2015
# PROFESSOR: HEDIBERT FREITAS LOPES
#
###########################################################################

###########################################################################
# Problema 1
###########################################################################
library(MASS)
n = 200
q = 2
A = matrix(c(0.6,0.1,0.8,0.8),q,q)
mu = c(0,0)
S = diag(100,q)

# a) 
set.seed(1234)
e = mvrnorm(n,mu,S)
y = matrix(0,n,q)
for (t in 2:n)
  y[t,] = A%*%y[t-1,]+e[t,]

# In my simualtion y(200) is equal to (-152,-101)

pdf(file="lista41.pdf",width=10,height=12)
par(mfrow=c(2,2))
plot(e,xlab="error (1st equation)",ylab="error (2nd equation)")
title("Uncorrelated errors")
plot(y,xlab="1st time series",ylab="2nd time series")
title("Bivariate VAR(1) data")
plot(y[,1],main="1st time series",ylab="",xlab="time",type="o")
plot(y[,2],main="2nd time series",ylab="",xlab="time",type="o")
dev.off()

# b)
hmax=20
As=array(0,c(hmax,q,q))
V=array(0,c(hmax,q,q))
yf=matrix(0,hmax,q)
As[1,,]=A
V[1,,] = S
yf[1,] = A%*%y[n,]
for (h in 2:hmax){
  As[h,,]=As[h-1,,]%*%A
  V[h,,] = V[h-1,,] + As[h-1,,]%*%S%*%t(As[h-1,,])
  yf[h,] = As[h,,]%*%y[n,]  
}  

pdf(file="lista42.pdf",width=10,height=8)
par(mfrow=c(1,2))
for (i in 1:q){
  L = yf[,i]-2*sqrt(V[,i,i])
  U = yf[,i]+2*sqrt(V[,i,i])
  ts.plot(cbind(L,yf[,i],U),col=c(2,1,2),lty=c(2,1,2),type="o")
}
dev.off()

# c)
# Como vimos na aula, uma forma alternativa de se verificar que o modelo e' estacionario
# e' obter os autovalores da matrix A.
# Os autovalores de A sao 1.0 e 0.4, ou seja, um deles esta' sobre o circulo unitario e, consequentemente,
# o processo e' nao-estacionario.
eigen(A)
$values
[1] 1.0 0.4

pdf(file="lista46.pdf",width=10,height=12)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
    plot(As[,i,j],type="o",pch=16,xlab="h",ylab="",main=paste("Phi(",i,",",j,")",sep=""))
    abline(h=0,col=2)
}
dev.off()




###########################################################################
# Problema 3
###########################################################################
library(MASS)

q=3
nu  = c(2,1,0)
A1 = matrix(c(0.7,0,0.9,0.1,0.4,0,0,0.1,0.8),q,q)
A2 = matrix(c(-0.2,0,0,0,0.1,0,0,0.1,0),q,q)
Sigma = matrix(c(0.26,0.03,0,0.03,0.09,0,0,0,0.81),q,q)

A = cbind(A1,A2)
I = diag(1,q)
Z = matrix(0,q,q)
A = rbind(A,cbind(I,Z))

mu = solve(I+A1+A2)%*%nu

h = 50
Phis = array(0,c(h,q,q))
Phis[1,,] = A1
Phis[2,,] = Phis[1,,]%*%A1+A2
for (i in 3:h)
Phis[i,,] = Phis[i-1,,]%*%A1+Phis[i-2,,]%*%A2

pdf(file="lista44.pdf",width=10,height=12)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
    plot(Phis[,i,j],type="o",pch=16,xlab="h",ylab="",main=paste("Phi(",i,",",j,")",sep=""))
    abline(h=0,col=2)
}
dev.off()


###########################################################################
# Problema 4
###########################################################################
library(MASS)
n    = 1000
q    = 2
A1 = matrix(c(0.5,0.4,0.1,0.5),q,q)
A2 = matrix(c(0.0,0.25,0.0,0.0),q,q)
mu = c(0,0)
S = diag(100,q)

set.seed(1234)
e = mvrnorm(n,mu,S)
y = matrix(0,n,q)
for (t in 3:n)
  y[t,] = A1%*%y[t-1,]+A2%*%y[t-2,]+e[t,]


pdf(file="lista43.pdf",width=12,height=8)
par(mfrow=c(2,1))
plot(y[,1],main="1st time series",ylab="",xlab="time",type="l")
plot(y[,2],main="2nd time series",ylab="",xlab="time",type="l")
dev.off()


h = 50
Phis = array(0,c(h,q,q))
Phis[1,,] = A1
Phis[2,,] = Phis[1,,]%*%A1+A2
for (i in 3:h)
Phis[i,,] = Phis[i-1,,]%*%A1+Phis[i-2,,]%*%A2

pdf(file="lista45.pdf",width=10,height=12)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
    plot(Phis[,i,j],type="o",pch=16,xlab="h",ylab="",main=paste("Phi(",i,",",j,")",sep=""))
    abline(h=0,col=2)
}
dev.off()