{

#test=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)
#test=test[,2:13]
#you have to string it to create a long vector
#x=array(t(test))

#sunspot data..
library(ts)
data(sunspot)
x=sunspot.year
N=length(x)

xeval=seq(min(x)-0.25*sd(x),max(x)+0.25*sd(x),length=100)
nevals=length(xeval)

library(sm)
zz=sm.density(x,eval.points=xeval,display="none")
xdensityorig=zz$estimate

#Fitting AR model
zz=ar(x,order.max=100,aic=TRUE)

#zz$order (gives you the best order)
#zz$aic   (gives the AIC values of models of different orders)
#zz$ar (gives the parameters of the best model based on AIC)
#zz$var.pred	(gives the error variance)

#Fitting an ARMA model:
#in R
#zz=arima0(x,order=c(1,0,1))	ARMA(1,1)
#zz=arima0(x,order=c(2,0,1))	ARMA(2,1)
#zz
#gives the all the model details
#also zz$aic gives the AIC

#in S
#zz=arima.mle(x,list(order=c(1,0,1)))	ARMA(1,1)
#etc.

#You will repeat this for different orders and pick the ARMA model
#that gives the minimum AIC
#------------------------

#Simulating from the fitted model:
#----------------------------------

nsim=100	#generate 100 simulates each of length N

simmean=1:nsim
simvar=1:nsim
simskew=1:nsim
simlag1=1:nsim
simpdf=matrix(0,nrow=nsim,ncol=nevals)

for(isim in 1:nsim){

zsim=arima.sim(n=N,list(ar=c(zz$ar)), sd=sqrt(zz$var.pred))+zz$x.mean

simmean[isim]=mean(zsim)
simvar[isim]=var(zsim)
simskew[isim]=skew(zsim)
simlag1[isim]=cor(zsim[1:(N-1)], zsim[2:N])
simpdf[isim,]=sm.density(zsim,eval.points=xeval,display="none")$estimate
}

par(mfrow=c(1,4))

zz=boxplot(simmean, plot=F)
z1=bxp(zz, ylim=range(simmean,mean(x)), xlab="",ylab="Mean")
points(z1,mean(x),pch=16,col="red")

zz=boxplot(simvar, plot=F)
z1=bxp(zz, ylim=range(simvar,var(x)), xlab="",ylab="Var")
points(z1,var(x),pch=16,col="red")

zz=boxplot(simskew, plot=F)
z1=bxp(zz, ylim=range(skew(x),simskew),xlab="",ylab="Skew")
points(z1,skew(x),pch=16,col="red")

zz=boxplot(simlag1, plot=F)
z1=bxp(zz, ylim=range(cor(x[1:(N-1)],x[2:N]),simlag1),xlab="",ylab="Skew")
points(z1,cor(x[1:(N-1)],x[2:N]),pch=16,col="red")

#boxplots of the PDF....
par(mfrow=c(1,1))
zz=boxplot(split(simpdf,col(simpdf)),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(simpdf,xdensityorig),xlab="",ylab="",cex=1.25)

z2=1:6
n1=1:6
z2[1]=z1[1]
z2[2]=z1[20]
z2[3]=z1[40]
z2[4]=z1[60]
z2[5]=z1[80]
z2[6]=z1[100]

n1[1]=xeval[1]
n1[2]=xeval[20]
n1[3]=xeval[40]
n1[4]=xeval[60]
n1[5]=xeval[80]
n1[6]=xeval[100]
n1=round(n1,dig=0)
n1=as.character(n1)

axis(1,at=z2,labels=n1,cex=1.00)
lines(z1,xdensityorig,lty=2,lwd=2,col="red")

}