##
    ## Fit a Seasonal Multivariate Autoregressive
    ##


    ## need to create 12 models one for each month.  Each one will be MAR(1)
    ## Z(t) = A * Z(t-1) + B Eps(t).

## Suppose the monthly streamflow matrix at Lees Ferry and Green River are
## in lf and gr, respectively.

## Lees Ferry   2006-1906
lf = matrix(scan("LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms

lfmean = colMeans(lf)
lfsd = apply(lf,2,sd)

lfscale = t( (t(lf) - lfmean)/lfsd )
#lfscale = scale(lf)

## Green River Ut
gr = matrix(scan("GreenRiver-UT-monflows-1906-2006.txt"), ncol=13, byrow=T)
gr = gr[,2:13]
gr = gr*0.0004690502     #convert to cms

grmean = colMeans(gr)
grsd = apply(gr,2,sd)

grscale = t( (t(gr) - grmean)/grsd )
#grscale = scale(gr)

    ## Bind the data into a 2 x 101(NumYears) X 12(months) matrix
library(abind)

    z = abind(lfscale, grscale, along=0)       # binds matrices.

nyrs = length(z[1,,1])   # number of years
nyrs1=nyrs-1

## Estimate the paramter matrices As and Bs for each month
As = array(list(), c(12))
Bs = array(list(), c(12))
for(imon in 1:12){
imon1=imon-1
if(imon == 1)imon1=12

if(imon == 1){

zt = t(z[,2:nyrs,imon])
zt1 = t(z[,1:nyrs1,imon1]) }
else
{
zt = t(z[,1:nyrs,imon])
zt1 = t(z[,1:nyrs,imon1]) }

M0 = var(zt)
M0.prev = var(zt1)
M1 = var(zt,zt1)

 ## solve for A
        ## A = M_1 * M_0^-1
        A = M1 %*% solve(M0.prev)
        As[imon] = list(A)

        ## Solve for B
        ## BB+ = D = M0 - M1 M0^-1 M1+
        D = M0 - M1 %*% solve(M0.prev) %*% t(M1)

        zz = svd(D)
        E = zz$u
        lambda = zz$d

        ## Check that I'm getting what I expected.  Yup this was zero
        ## D - E %*% diag(lambda) %*% t(E)

        B = E %*% sqrt(diag(lambda))
        Bs[imon] = list(B)


}

########### simulate..
nsims=250   # number of simulation..
sims = array(0,dim=c(nsims,dim(z)))

### sims is a four dimensional array..
### first dimension is the simulation number
### second dimension is the location 1, for lees ferry and 2 for green river
### third dimension is rows - i.e., the years
### fourth is the columns - i.e, months
for(i in 1:nsims){

this.sim = sims[i,,,]

## Start with a previ Dec value..
sample.num = round(runif(1,1,nyrs))
Z.prev = as.matrix(z[,sample.num,12])

for(j in 1:nyrs){
 
    for(imon in 1:12){

        epsilon = as.matrix(rnorm(2,0,1))     #normal random 

        Z.t = As[[imon]] %*% Z.prev + Bs[[imon]] %*% epsilon

## Assign the predicted value.

	this.sim[,j,imon] = Z.t

## set the previous value and continue.
	Z.prev = Z.t
			}
		}

## Put back the mean and std dev..  and 
## Store the simulation back in the big sim array

this.sim[1,,] = t((t(this.sim[1,,]) * lfsd) + lfmean)

this.sim[2,,] = t((t(this.sim[2,,]) * grsd) + grmean)

sims[i,,,] = this.sim

}

#### You can grab each simulation for each location separately and compute
### the stats and do the boxplots. You can use the codes you have from 
### HW 2.