# Simulate and review a Bivariate AR(2)

phione<-matrix(nrow=2,byrow=T,c(.5,.9,0,.3))
phitwo<-matrix(nrow=2,byrow=T,c(0,0,0,.4))

phi<-array(0,dim=c(2,2,2))
phi[,,1]<-phione
phi[,,2]<-phitwo

#compute infinite MA infinity coefficients

	identity<-matrix(nrow=2,byrow=T,c(1,0,0,1))
	m<-40				#truncation point
	psi<-array(0,dim=c(2,2,m))

	psi[,,1]<-identity
	psi[,,2]<-phi[,,1]

	for (i in 3:m) psi[,,i]<-phi[,,1]%*%psi[,,i-1] + phi[,,2]%*%psi[,,i-2]
	
#Simulate process

	n<-100
	npm<-n+m
	z<-matrix(0,nrow=n,ncol=2)
	a<-matrix(rnorm(npm*2,0,1),nrow=2, ncol=npm)
	temp<-matrix(0,nrow=2,ncol=1)
	for(i in 1:n) {
		for (j in 1:m) z[i,]<-z[i,]+psi[,,j]%*%a[,i+j]
	}

# Preliminary Plots

graphsheet()
	par(mfrow=c(3,2),omi=c(.5,1,.5,1))

	tsplot(psi[1,1,],psi[1,2,],psi[2,1,],psi[2,2,],type="l")
	title(main="MA Infinity Coef")
	tsplot(z,type="l",lty=c(1,1),col=c(1,4))
	title(main="Time Series Plot of Each Series")
	acf(z)
	acf(z,type="partial")

# Fit a Bivariate AR

	zarfit<-ar(z,aic=T,order.max=2,method="burg")
	zarfit$order
	zarfit$ar
	zarfit$var.pred
	tsplot(zarfit$resid,type="l",lty=c(1,1),col=c(1,4))
	title(main="Residuals After Bivariate AR Fit")
	acf(zarfit$resid)

# Plot Estimated Spectral Density

graphsheet()
	par(mfrow=c(2,2),omi=c(.5,1,.5,1))

	bivspec<-spec.ar(zarfit,plot=T)  #(z,plot=T,method="AR")	
	plot(bivspec$freq,bivspec$coh,type="l")
	title(main="Est. Coherencey Spectrum")
	plot(bivspec$freq,bivspec$phase,type="l")
	title(main="Est. Phase Spectrum")

	bivspec<-spectrum(z,plot=F,method="pgram",spans=10)	
	matplot(bivspec$freq,cbind(bivspec$spec+1,bivspec$coh,bivspec$phase),type="l")
	title(main="Smoothed Sample Bivariate Spectrum")

# Plot True Spectral Density

zartrue<-zarfit
zartrue$ar[1,,]<- -1*phi[,,1]
zartrue$ar[2,,]<- -1* phi[,,2]
zartrue$order<-2


graphsheet()
	par(mfrow=c(2,2),omi=c(.5,1,.5,1))

	bivspec<-spec.ar(zartrue,plot=T)  
	plot(bivspec$freq,bivspec$coh,type="l")
	title(main="Squared Coherencey Spectrum")
	plot(bivspec$freq,bivspec$phase,type="l")
	title(main="Phase Spectrum")