source("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/ssa-b.txt")

#### This code does the SSA and produces the RCs, plots the fraction
### Eigen values (similar to the Eigen spectrum of PCa)
### the first part of problem 7

## For forecasting, fit AR models to the RCs
## Suppose RC is the reconstructed component through time t and you want 
## a one step ahead forecast
#
#rc.ar = ar(RC)
#xp = predict(rc.ar,n.ahead=1)
#xp$pred gives the prediction
#xp$se gives the prediction standard error

#Repeat this for all the RCs. As before, fit a Normal distribution to the
#noise RCs

## The mean forecast is the sum of all the component forecasts
### The standard error of this is = sqrt(serc1^2 + serc2^2 + ... + serck^2)
## Assuming normal distribution you can compute the threshold probabilities for
## RPSS


source("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session4/ssa-b.txt")
lf = matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session4/LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
year = lf[,1]
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms
lfann = apply(lf,1,sum)

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

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


### 101 years of monthly data - M should be roughtly N/5 = ~ 20years
#M = 20*12 	#months
#ssaout = ssab(lfscale,M)


### SSA on Annual total volume
M = 20
ssaout = ssab(lfann,M)

#ssaout$lambdas = Fraction Eigen Value
## ssaout$Rpc = reconstructed components


## Add the first K reconstructed components
K = 7		# just as an example

Recon = apply(ssaout$Rpc[,1:K],1,sum)

plot(year, lfann, type="l", xlab="Year", ylab="Lees Ferry Annual Flow", ylim=range(lfann,Recon))
lines(year, Recon, col="red", lwd=2)

## as a check if you sum *all* the reconstructed components you should get the original
### data back

xrecover = apply(ssaout$Rpc,1,sum)

## this should be equal to lfann


### Get the seasonal NINO34 data
(Dec 1856 - Jan 1857) to (Sep 2009 - Nov 2009)
153 years and 4 seasons each year.
nino34ses=scan("http://iridl.ldeo.columbia.edu/expert/SOURCES/.Indices/.nino/.EXTENDED/.NINO34/T/%28Dec%201856%29%28Nov%202009%29RANGE/T/3/boxAverage/data.ch")
M = round (10 * 4)		#10-years

ssaout = ssab(nino34ses,M)

K = 10
Recon = apply(ssaout$Rpc[,1:K],1,sum)
plot(ts(nino34ses,start=1857,frequency=4,end=2009), xlab="Year", ylab="Seasonal NINO34",
ylim=range(nino34ses,Recon))
lines(ts(Recon,start=1857,frequency=4,end=2009),lwd=2, col="red")