{

#This code fits a disaggregation model to the monthly and annual time series
# then simulates from the fitted model..

#Y is the matrix of n x m (n is the number of years and m is the number of
# seasons, for monthly data m = 12)

#X is a vector of length n (which are the annual values)

data=matrix(scan("Leesferry-mon-data.txt"),ncol=13,byrow=T)
data=data[,2:13]/10000                   #remove year and divide
#y=data[,1:2]
x=rowSums(y)
N=length(x)
M=length(data[1,])

Sxx=var(x)	#gives one value 1x1
Syy=var(y)	#gives a m x m matrix
Syx=var(y,x)	#gives a m x 1 vector
Sxy=var(x,y)	#gives a 1 x m vector

A=Syx/Sxx	#gives a m x 1 vector

D=Syy - (Syx/Sxx) %*% Sxy	#gives a m x m matrix
zz=svd(D)
B=zz$u %*% diag(sqrt(zz$d)) 	#gives a m x m matrix

C=rep(1,M)

#check
zz=C %*% B	#(this should be a matrix of 1 x m and should be equal to 0)

zz=C %*% A      #(this should be equal to 1)

#now for simulation..
#let us suppose that you have a model (e.g. AR(1)) on the annual
#values to generate annual values
#let us suppose that you have the simulated annual values of
#length n years in xsim (n x 1)
W=matrix(rnorm(1,0,0),nrow=12,ncol=N)
annmodel=ar(x-mean(x))
xsim=arima.sim(n=N,list(ar=c(annmodel$ar)), sd=sqrt(annmodel$var.pred))+mean(x)


xsim=A %*% t(xsim) + B %*% W

xsim=t(xsim)	#this gives a n x m matrix of simulated values from
		#the disaggregagted model. One can now compute the
		#the statistics of the various seasons etc.


}