prettyplot<-function(nrow,ncol){
	par(mfrow=c(nrow,ncol),omi=c(.5,1,.5,1),lwd=2)
}

prettyplot(2,1)

# AR spectral density

armodel<-list(ar=c(0,0,-.8),order=3,var.pred=1)
arspec<-spec.ar(armodel,n.freq=100,plot=T)

# MA spectral density

spec.ma<-function(ma,order,var.pred,n.freq){
	maspec<-spec.ar(list(ar=ma,order=order,var.pred=1),n.freq=n.freq,plot=F)
	temp<-10^(maspec$spec/10)  # reverse the transformations see spec.ar
	temp2<-var.pred/temp  # properly include the variance 
	maspec$spec<-10*log10(temp2)  # perform the same transf. as in spec.ar
	plot(maspec$freq,maspec$spec,type="l")
	return(maspec)
}

maspec<-spec.ma(c(.9,0,-.5),3,1,100)

# ARMA spectral density
spec.arma<-function(ar,ma,order,var.pred,n.freq){
	print(list(ar=ar,ma=ma,order=order,var.pred=var.pred,n.freq=n.freq))
	prettyplot(3,1)
	#compute AR component
	arspec<-spec.ar(list(ar=ar,order=order[1],var.pred=1),n.freq=n.freq,plot=T)
	temp<-10^(arspec$spec/10) # reverse the transformation
	
	#compute MA component
	maspec<-spec.ar(list(ar=ma,order=order[2],var.pred=1),n.freq=n.freq,plot=F)
	temp1<-10^(maspec$spec/10)  # reverse the transformations see spec.ar
	temp2<-var.pred/temp1  # properly include the variance 
	maspec$spec<-10*log10(temp2)  # perform the same transf. as in spec.ar
	plot(maspec$freq,maspec$spec,type="l")
	title("Spectral Density of MA Component")

	
	#put it all together
	temp2<-var.pred*temp/temp1  
	spec<-10*log10(temp2)  # perform the same transf. as in spec.ar
	armaspec<-list(spec=spec,freq=arspec$freq,order=order,var.pred=var.pred)
	plot(armaspec$freq,armaspec$spec,type="l")
	title("ARMA Spectral Density")
	return(armaspec)
}

armaspec<-spec.arma(c(.9,0,-5),c(0,0,-.2),c(3,3),1,100)