#simulate from a PDF, compute the stats for each Monte Carlo simulation
# and boxplot them along with those of the historical data
# Example is with lognormal. You can see the commands to do
### the same with any other PDF

library(MASS)

data=matrix(scan("aismr.txt"),ncol=13,byrow=T)
xdata=data[,2]	#Jan Rainfall

#suppose the historical data is in the vector 'xdata'
N=length(xdata)		#number of data points.

#Monte Carlo...
nsim=1000	#number of simulations, each of length N

#points at which to estimate the PDF from the simulations and the observed

xeval=seq(min(xdata)-sd(xdata),max(xdata)+sd(xdata),length=100)
neval=length(xeval)
xsimpdf=matrix(0,nrow=nsim,ncol=neval)

#initialize vectors for the statistics
meansim=1:nsim
mediansim=1:nsim
sdsim=1:nsim
skewsim=1:nsim
iqrsim=1:nsim

#### Select a model (i.e., PDF) to fit to the data and 
# simulate from. Typically, you will select this based on your
# diagnostics and goodness of fit exercise - using visual and KS tests

### If you decide on Gamma or Weibull you have to fit it to the data
# to obtain the paramteres..

zgamma=fitdistr(xdata,"gamma")
zweibull=fitdistr(xdata,"weibull")

###########  Estimate the PDF of the selected model at athe evaluation points 
# based on the data.

## kernel PDf
xdensitykernel = sm.density(xdata,eval.points=xeval,display="none")$estimate

#lognormal
#xdensityorig=dlnorm(xeval,meanlog=mean(log(xdata)), sdlog=sd(log(xdata)))

#normal
#xdensityorig=dnorm(xeval,meanlog=mean(xdata), sdlog=sd(xdata))

#exponential
#xdensityorig=dexp(xeval,rate=1/mean(xdata))

# Gamma
xdensityorig=dgamma(xeval,shape=zgamma$estimate[1],scale=1/zgamma$estimate[2])

# Weibull
#xdensityorig=dweibull(xeval,shape=zweibull$estimate[1],scale=zweibull$estimate[2])

######## Simulation ##########

for(i in 1:nsim){

# Simulate from the desired model (i.e., PDF) 
# each simulation should of the same length as the original data
# used to fit the model, N = length(xdata)

#simulate from Lognormal
#xsim=rlnorm(N,meanlog=mean(log(xdata)), sdlog=sd(log(xdata)))

#simulate from Normal
#xsim = rnorm(N, mean=mean(xdata), sd=sd(xdata))

#Simulate from Exponential..
#xsim = rexp(N, rate=1/mean(xdata))


#simulate from Gamma..
xsim=rgamma(N,shape=zgamma$estimate[1],scale=1/zgamma$estimate[2])

#simulate from Weibull..
#xsim=rweibull(N,shape=zweibull$estimate[1],scale=zweibull$estimate[2])


# compute the statistics from the simulation
meansim[i]=mean(xsim)
sdsim[i]=sd(xsim)
#sdsim[i]=sd(log(xsim))
skewsim[i]=skew(xsim)
iqrsim[i]=diff(quantile(xsim,c(0.25,0.75)))
mediansim[i]=quantile(xsim,c(0.5))

#estimate the PDF at the evaluation points based on the simulated data
#replace with the appropriate distribution..

#Log Normal..
xsimpdf[i,]=dlnorm(xeval,meanlog=mean(log(xsim)), sdlog=sd(log(xsim)))

#Exponential...
#xsimpdf[i,]=dexp(xeval,rate=1/mean(xsim))

#Normal 
#xsimpdf[i,]=dnorm(xeval,meanlog=mean(xsim), sd=sd(xsim))

#Gamma..
#zz=fitdistr(xsim,"gamma")
#xsimpdf[i,]=dgamma(xeval,shape=zz$estimate[1],scale=1/zz$estimate[2])

#Weibull..
#zz=fitdistr(xsim,"weibull")
#xsimpdf[i,]=dweibull(xeval,shape=zz$estimate[1],scale=zz$estimate[2])
}


#boxplot the stats from the simulations and overlay the corresponding value
# from the original data as a point.

par(mfrow=c(1,5))

zz=boxplot(meansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,mean(xdata),col="red",cex=1.25,pch=16)
title(main="Mean")

zz=boxplot(mediansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,quantile(xdata,c(0.5)),col="red",cex=1.25, pch=16)
title(main="Median")

zz=boxplot(sdsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
#points(z1,sd(log(xdata)),col="red",cex=1.25)
points(z1,sd(xdata),col="red",cex=1.25,pch=16)
title(main="SD")

zz=boxplot(skewsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,skew(xdata),col="red",cex=1.25,pch=16)
title(main="Skew")

zz=boxplot(iqrsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,diff(quantile(xdata,c(0.25,0.75))),col="red",cex=1.25,pch=16)
title(main="IQR")

par(ask=TRUE)

#boxplot the PDFs from the simulations along with that of the original data

par(mfrow=c(1,1))
xs=1:neval

# For R version  2.4.1 
zz=boxplot(split(t(xsimpdf),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(xsimpdf,xdensityorig),xlab="",ylab="",cex=1.25)

npts=10         #number of points to plot on the x-axis..
n2=round(neval*(1:npts)/npts)

z2=z1[n2]
n1=xeval[n2]


n1=round(n1,dig=2)
n1=as.character(n1)

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

lines(z1,xdensitykernel,col="blue",lwd=2)

title(main="PDFs from the simulations and the historical data")