# 2014-03-31
# start with this script (a)
# show script (b) version with matrix notation
# example of covariance effecting the variance
# in the averaging
# clear memory
rm(list=ls())
# reset the plot window to full sized plot
par(mfrow=c(1,1))
source("scripts/loadano2.R")



month<-"Jan"
# station 1
x1<-loadano2(station="USW00014735",state="NY",month=month,start=1950,end=2012)
# Station 2
x2<-loadano2(station="USW00014750",state="NY",month=month,start=1950,end=2012)
# Station 3
x3<-loadano2(station="USW00014771",state="NY",month=month,start=1950,end=2012)
# Station 4
x4<-loadano2(station="USW00094728",state="NY",month=month,start=1950,end=2012)


n<-length(x1$ano)

xm1<-mean(x1$ano)
xm2<-mean(x2$ano)
xm3<-mean(x3$ano)
xm4<-mean(x4$ano)

x1$ano<-x1$ano-xm1
x2$ano<-x2$ano-xm2
x3$ano<-x3$ano-xm3
x4$ano<-x4$ano-xm4


z<-(x1$ano+x2$ano+x3$ano+x4$ano)

szz<-round(var(z),4)

print("Variance of the variable z=x1+x2+x3+x4")
print(paste("Var(z)    : ",szz,sep=''))



sxx1<-round(var(x1$ano),4)
sxx2<-round(var(x2$ano),4)
sxx3<-round(var(x3$ano),4)
sxx4<-round(var(x4$ano),4)

print ("Variance of the anomalies")
print (paste("station 1 : ",sxx1,sep=''))
print (paste("station 2 : ",sxx2,sep=''))
print (paste("station 3 : ",sxx3,sep=''))
print (paste("station 4 : ",sxx4,sep=''))
print (paste("---------------------------"))
ssum<-sxx1+sxx2+sxx3+sxx4
print (paste(" summed up:",ssum,sep=''))





       
# adding the covariances
c12<-round(cov(x1$ano,x2$ano),4)
c13<-round(cov(x1$ano,x3$ano),4)
c14<-round(cov(x1$ano,x4$ano),4)
c23<-round(cov(x2$ano,x3$ano),4)
c24<-round(cov(x2$ano,x4$ano),4)
c34<-round(cov(x3$ano,x4$ano),4)

print ("Covariance terms for the anomalies")
print (paste("station 1,2 : ",c12,sep=''))
print (paste("station 1,3 : ",c13,sep=''))
print (paste("station 1,4 : ",c14,sep=''))
print (paste("station 2,3 : ",c23,sep=''))
print (paste("station 2,4 : ",c24,sep=''))
print (paste("station 3,4 : ",c34,sep=''))
print (paste("---------------------------"))
csum<-c12+c13+c14+c23+c24+c34
print (paste(" summed up:",csum,sep=''))

tsum<-ssum+2*csum
print (paste("sum of variance and 2*covariance terms : ",tsum))