# PROGRAM 10
#
# K. Ensor (ensor@rice.edu)
# Rice Univeristy
#

#DESCRIPTIVE PLOTS - rewriting this to include the histogram
# 						 and the default acf/pacf

kdescripc<-function(x,xtitle)
{
split.screen(c(2,1)) #split screen into two parts
split.screen(c(1,3),screen=2) #split the bottom graph into three parts

screen(1)
tsplot(x)
title(main=xtitle)

screen(3)
hist(x)

screen(4)
acf(x)

screen(5)
acf(x,type="partial")

close.screen(all=T)
}

# K. Ensor
# Step 2 in ARIMA diagnostics
#
karima.diag<-function(stresid,varest,k){
	n=length(stresid)
	par(mfrow=c(2,1))
	qqnorm(stresid)
	hist(stresid)
	par(mfrow=c(1,1))
	
	loglik<-log(varest)
	aic<-loglik+ 2*k/n
	aicc<-loglik+(n+k)/(n-k-2)
	sic<-loglik+ k*log(n)/n
	return(aic,aicc,sic)		
}

# Uses the seasonally adjusted US GNP referred to as: ss.gnp

kdescripc(ss.gnp,"Seasonally Adjusted U.S. GNP from 1947(1) to 1991(1)")

# Trend hides everything
# Let's remove the trend

x<-diff(ss.gnp,1)
kdescripc(x,"First Difference of SA GNP")

# We know see that the variance increase with time,
# so let's go back and take the log of the data and
# then difference

ss.lgnp<-log(ss.gnp)
kdescripc(ss.lgnp,"Log of SA GNP")

x<-diff(ss.lgnp,1)
kdescripc(x,"First Difference of LSA GNP")

# Two different models MA(2) and AR(3) - why?

gnpma<-arima.mle(x,model=list(order=c(0,0,2)),xreg=1)
gnpmad<-arima.diag(gnpma)
gnpmadk<-karima.diag(gnpmad$std.resid,gnpma$sigma2,2)

gnpar<-arima.mle(x,model=list(order=c(3,0,0)),xreg=1)
gnpard<-arima.diag(gnpar)
gnpardk<-karima.diag(gnpard$std.resid,gnpar$sigma2,3)

################## STOP HERE
#########################################