{

source("skew.s")
# This code is to model a Seasonal AR(1) or Thomas-Fiering Model
# for monthly streamflow time series.

# you should HAVE THE DATE IN THIS MATRIX WYmonflow before you 
# run this code..

#the data is in WYmonflow - that contains number of years as
# rows and 12 columns (Jan through Dec)


test=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)
WYmonflow=test[,2:13]
WYmonflow=WYmonflow/10000


nsim=250	# number of simulations 

nsim1=nsim+1

#initialize the matrices that store the statistics
# the first row contains the statistic of the historical
#data. 

armean=matrix(0,nsim1,12)	
arstdev=matrix(0,nsim1,12)
arcor=matrix(0,nsim1,12)
arskw=matrix(0,nsim1,12)
armax=matrix(0,nsim1,12)
armin=matrix(0,nsim1,12)

nyrs=length(WYmonflow[,1])	#number of years of data for each month
nyrs1=nyrs-1


# get the co-efficients of the Thomas Fiering Model 
coef1=1:12
coef2=1:12
for(j in 1:12){
j1=j-1
if(j == 1)j1=12 
if(j == 1)coef1[j]=(cor(WYmonflow[2:nyrs,j], WYmonflow[1:nyrs1,j1])) 
if(j > 1)coef1[j]=cor(WYmonflow[1:nyrs,j], WYmonflow[1:nyrs,j1])
coef2[j]=sqrt((var(WYmonflow[,j])) * (1. - coef1[j]*coef1[j]))
coef1[j]=coef1[j]*sqrt(var(WYmonflow[,j])/var(WYmonflow[,j1]))
}


for(k in 1:nsim){
nmons=nyrs*12	#number of values to be generated
		#nyrs is the number of years and 12 is the number of
		#months.

xsim=1:nmons

i=round(runif(1,1,nyrs))
xsim[1]=WYmonflow[i,1]
xprev=xsim[1]

for(i in 2:nmons){
j=i %% 12
if(j == 0)j=12
j1=j-1

if(j == 1)j1=12

x1=xprev-mean(WYmonflow[,j1])
x2=coef2[j]*rnorm(1,0,1)

xsim[i]=mean(WYmonflow[,j]) + x1*coef1[j] + x2

xprev=xsim[i]}

#put the array inot a matrix 
simdismon=matrix(xsim,nrow=12)
simdismon=t(simdismon)

#***calculate basic statistics****************

k1=k+1
for(j in 1:12){

                armean[k1,j]=mean(simdismon[,j])
                armax[k1,j]=max(simdismon[,j])
                armin[k1,j]=min(simdismon[,j])
                arstdev[k1,j]=sd(simdismon[,j])
		arskw[k1,j]=skew(simdismon[,j])

#                arskw[k1,j]=sum((simdismon[,j]-armean[k1,j])^3)
#                arskw[k1,j]=arskw[k1,j]/21
#                arskw[k1,j]=arskw[k1,j]/arstdev[k1,j]^3
        }
for(j in 2:12) {
j1=j-1
                arcor[k1,j]=cor(simdismon[,j],simdismon[,j1])
        }
arcor[k1,1]=cor(simdismon[1:nyrs1,12],simdismon[2:nyrs,1])

}

#*************************
#Basic Data

obsmean=1:12
obsstdev=1:12
obscor=1:12
obsskw=1:12
obsmax=1:12
obsmin=1:12

for(i in 1:12){

        obsmax[i]=max(WYmonflow[,i])
        obsmin[i]=min(WYmonflow[,i])
        obsmean[i]=mean(WYmonflow[,i])
        obsstdev[i]=sd(WYmonflow[,i])
	obsskw[i]=skew(WYmonflow[,i])

}
for(i in 2:12){
i1=i-1
        obscor[i]=cor(WYmonflow[,i], WYmonflow[,i1])
}
obscor[1]= cor(WYmonflow[1:nyrs1,12], WYmonflow[2:nyrs,1])

        armean[1,1:12]=obsmean[1:12]
        arstdev[1,1:12]=obsstdev[1:12]
        arskw[1,1:12]=obsskw[1:12]
        arcor[1,1:12]=obscor[1:12]
        armax[1,1:12]=obsmax[1:12]
        armin[1,1:12]=obsmin[1:12]

#write the stats in separate files..

write(t(arcor),file="arcorrs1",ncol=12)
write(t(armean),file="armeans",ncol=12)
write(t(arstdev),file="arstdevs",ncol=12)
write(t(arskw),file="arskews",ncol=12)
write(t(armax),file="armax",ncol=12)
write(t(armin),file="armin",ncol=12)

}