#Run from command line using R CMD BATCH ar1_example.r 
#Run from inside R using source("ar1_example.r",echo=T)

num <- 100
desired_width <- 800
desired_height <- 800
desired_res <- 120

legends <- c("phi=0.0", "phi=0.2", "phi=0.9","Trend + 'phi=0.2'")
pt_two_col <- "forestgreen"
colors <- c("black", pt_two_col, "red", pt_two_col)
width <- 3
widths <- c(width,width,width,width)
ltys <- c(2,1,6,3)
actualvalues <- c(0.0,0.2,0.9,0.2)
white   <- arima.sim(n = num, list(ar = 0.0), innov=rnorm(num))
pt_two  <- arima.sim(n = num, list(ar = 0.2), innov=rnorm(num))
pt_nine <- arima.sim(n = num, list(ar = 0.9), innov=rnorm(num))
signal_plus_pt_two <- pt_two + 0.1*(1:num)
white   <- white   - mean(white)
pt_two  <- pt_two  - mean(pt_two)
pt_nine <- pt_nine - mean(pt_nine)
signal_plus_pt_two <- signal_plus_pt_two - mean(signal_plus_pt_two)
png(file="fig1.png",width=desired_width,height=desired_height,res=desired_res)
ts.plot(white, pt_two, pt_nine, signal_plus_pt_two, gpars=list(xlab="time", ylab="", main="AR(1) noise time series with different parameters 'phi'", col=colors, lty=ltys, lwd = widths))
legend("topleft", legend = legends, col=colors, lty=ltys, lwd=widths, cex=1)
dev.off

print(.Random.seed)

white_noise <- arima(white,order=c(1,0,0))
pt_two_noise <- arima(pt_two,order=c(1,0,0))
pt_nine_noise <- arima(pt_nine,order=c(1,0,0))
signal_plus_noise <- arima(signal_plus_pt_two,order=c(1,0,0))
print(white_noise)
print(pt_two_noise)
print(pt_nine_noise)
print(signal_plus_noise)

png(file="fig2.png",width=desired_width,height=desired_height,res=desired_res)
plot(acf(white,plot=FALSE), type="l", main="Autocorrelation function (ACF) of AR(1) time series", max.mfrow=1,col=colors[1], xlab="lag", ylab="", lty=ltys[1], lwd = widths[1])
lines(acf(pt_two,plot=FALSE)$acf[-1,1,1], lty=ltys[2], col=colors[2], lwd=widths[2])
lines(acf(pt_nine,plot=FALSE)$acf[-1,1,1], lty=ltys[3], col=colors[3], lwd=widths[3])
lines(acf(signal_plus_pt_two,plot=FALSE)$acf[-1,1,1], lty=ltys[4], col=colors[4], lwd=widths[4])
legend("topright", legend = legends, col=colors, lty=ltys, lwd=widths, cex=1)
dev.off

means <- c(white_noise$coef[1], pt_two_noise$coef[1], pt_nine_noise$coef[1], signal_plus_noise$coef[1])
names <- legends
standardErrors <- c(sqrt(diag(vcov(white_noise))) [1], sqrt(diag(vcov(pt_two_noise))) [1], sqrt(diag(vcov(pt_nine_noise))) [1], sqrt(diag(vcov(signal_plus_noise))) [1])

#Plot actual and estimated AR(1) parameters.
offset <- 0.2
width <- 1
separation <- 1.2
actualhalfwidth <- 0.02
lowest <- min(c(0,min(means - 2*standardErrors))) - offset/2
highest <- max(means + 2*standardErrors)
plotx <- offset + c(0,width,width,0,NA,separation,separation+width,separation+width,separation,NA,2*separation,2*separation+width,2*separation+width,2*separation,NA,3*separation,3*separation+width,3*separation+width,3*separation)
ploty1 <- c(means[1]+2*standardErrors[1],means[1]+2*standardErrors[1],means[1]-2*standardErrors[1],means[1]-2*standardErrors[1],NA,means[2]+2*standardErrors[2],means[2]+2*standardErrors[2],means[2]-2*standardErrors[2],means[2]-2*standardErrors[2],NA,means[3]+2*standardErrors[3],means[3]+2*standardErrors[3],means[3]-2*standardErrors[3],means[3]-2*standardErrors[3],NA,means[4]+2*standardErrors[4],means[4]+2*standardErrors[4],means[4]-2*standardErrors[4],means[4]-2*standardErrors[4])
ploty2 <- c(actualvalues[1]+actualhalfwidth,actualvalues[1]+actualhalfwidth,actualvalues[1]-actualhalfwidth,actualvalues[1]-actualhalfwidth,NA,actualvalues[2]+actualhalfwidth,actualvalues[2]+actualhalfwidth,actualvalues[2]-actualhalfwidth,actualvalues[2]-actualhalfwidth,NA,actualvalues[3]+actualhalfwidth,actualvalues[3]+actualhalfwidth,actualvalues[3]-actualhalfwidth,actualvalues[3]-actualhalfwidth,NA,actualvalues[4]+actualhalfwidth,actualvalues[4]+actualhalfwidth,actualvalues[4]-actualhalfwidth,actualvalues[4]-actualhalfwidth)
png(file="fig3.png",width=desired_width,height=desired_height,res=desired_res)
barCenters <- barplot(means, names.arg=names, cex.names=1.2, yaxt="n", border=NA, col=NA, las=1, ylim=c(lowest,highest))
axis(side = 2, at = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0), las=1)
title(main = "Actual vs estimated AR(1) parameters 'phi'")
polygon(plotx, ploty1, density = NA, border = "grey", col = "grey")
polygon(plotx, ploty2, density = NA, border = colors, col = colors)
dev.off