{

#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))

source("skew.s")

#sunspot data..
data(sunspot)
x=sunspot.year
N=length(x)
xmean=mean(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 Np-AR model

x=x-mean(x)

#based on MI we identify the lag to be

ilag=3

N1=N-ilag-1
xdata=matrix(0,ncol=ilag,nrow=N1)
ydata=1:N1

for(i in 1:N1){
i1=i
i2=i+ilag-1
i3=i2+1
xdata[i,]=x[i1:i2]
ydata[i]=x[i3]
}

K=round(sqrt(N1))
#K=K/2
W=1:K
W=1/W
W=W/sum(W)
W=cumsum(W)

nsim=50	#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)

zsim=1:N

for(isim in 1:nsim){

#first 1:lag points select at random

i=round(runif(1,1,N1))
xp=x[i:(i+ilag-1)]
zsim[1:ilag]=xp

for(j in (ilag+1):N){

xx=rbind(xp,xdata)
xdist=order(as.matrix(dist(xx))[1,2:(N1+1)])
xx=runif(1,0,1)
xy=c(xx,W)
xx=rank(xy)
i1=xdist[xx[1]]

zsim[j]=ydata[i1]

xp=zsim[(j-ilag+1):j]
}
zsim=zsim+xmean


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,xmean), xlab="",ylab="Mean")
points(z1,xmean,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")

dev.copy(postscript,file="npsim-stat-boxplots.ps")
dev.off()

#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")

dev.copy(postscript,file="npsim-pdf-boxplots.ps")
dev.off()
}