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")