Austria Warming 
 
 1 Global Warming
 from   "http://en.wikipedia.org/wiki/Global_warming" 
  Global warming is the rising average temperature of Earth's atmosphere and oceans since the late 19th century and its projected continuation. Since the early 20th century, Earth's average surface temperature has increased by about 0.8 °C (1.4 °F), with about two thirds of the increase occurring since 1980. Warming of the climate system is unequivocal, and scientists are more than 90% certain that most of it is caused by increasing concentrations of greenhouse gases produced by human activities such as deforestation and the burning of fossil fuels. ..Climate model projections .. indicate that during the 21st century the global surface temperature is likely to rise a further 1 to 6 °C . The ranges of these estimates arise from the use of models with differing sensitivity to greenhouse gas concentrations. Although our output of 29 gigatons of CO2 is tiny compared to the 750 gigatons moving through the carbon cycle each year, it adds up because the land and ocean cannot absorb all of the extra CO2. About 40% of this additional CO2 is absorbed. The rest remains in the atmosphere, and as a consequence, atmospheric CO2 is at its highest level in 15 to 20 million years (Tripati 2009). (A natural change of 100ppm normally takes 5,000 to 20,000 years. The recent increase of 100ppm has taken just 120 years). An increase in global temperature will cause sea levels to rise ..Warming is expected to be strongest in the Arctic and would be associated with continuing retreat of glaciers, permafrost and sea ice. …..  If global mean temperature increases to 4°C above pre-industrial levels, the limits for human adaptation are likely to be exceeded in many parts of the world, while the limits for adaptation for natural systems would largely be exceeded throughout the world. Hence, the ecosystem services upon which human livelihoods depend would not be preserved.
 2 Worldwide weather station locations 
 from http://cdiac.ornl.gov/epubs/ndp/ndp041/graphics/ndp041.temp.gif 
  There are at present about 40000 weather stations worldwide, concentrated in the U.S., Europe, and east Asia. There are few long established stations, and the more recent stations offer no evidence that the present increase in temperature is unusual compared to, say, 19th century values. There were perhaps only 100 weather stations in 1850. Berlin had one in 1701, Vienna had one in 1775, and Austria had 10 in 1850. We will look at what happened at these Austrian stations. We find data and plots from many long established European stations at: 
"http://www.zamg.ac.at/histalp/content/view/35/1/index.html" 
 3  Vienna and some long established German stations 
 picture( 'p/g2.png') 
picture( 'p/g3.png') 
picture( 'p/g4.png') 

    

After declining about a degree(C) in the 19th century, temperatures increased about two degrees over the last century, more sharply recently.
 4 Find historic Austrian temperature series 
 4 Find historic Austrian temperature series 
 We use temperature series going back to 1858, available in directory 
 C:/Users/John/Desktop/_ViennaData/Austrian Temperatures/data
  
dd <- 'C:/Users/John/Desktop/_ViennaData/Austrian Temperatures/data/' 
filelist <- list.files(dd) 
print(filelist) 
 [1] "boundary.txt"            "HISTALP-T01-GRA-hom.txt" "HISTALP-T01-INN-hom.txt" "HISTALP-T01-KLA-hom.txt" "HISTALP-T01-KRE-hom.txt"
 [6] "HISTALP-T01-LIN-hom.txt" "HISTALP-T01-OBG-hom.txt" "HISTALP-T01-SAL-hom.txt" "HISTALP-T01-SAN-hom.txt" "HISTALP-T01-VIA-hom.txt"
[11] "HISTALP-T01-WIE-hom.txt"

 5 Read data and select desired years  
 We make a data frame 10 x 154 to contain the 10 stations, 
 4 descriptive variables, and 150 yearly temperatures. 
  
tempD <- data.frame(matrix(0,10,154) ) 
names(tempD) <-c('loc', 'long', 'lat', 'alt', 1859:2008) 
# Run through all files but the first, which has Austrian boundary data. 
for ( ff in 2:length(filelist) ){ 
nf <- ff-1 
tf <-scan(paste(dd, filelist[ff], sep=''), what='', sep='') 
if (nf==1){ 
print(head(tf)) 
} 
 
# get station name, title, lat long altitude 
tempD[nf, 1] <- substr(tf[3], 1, 25) 
tempD[nf, 2] <- as.numeric(substr(tf[3], 36, 40)) 
tempD[nf, 3] <- as.numeric(substr(tf[3], 42, 46)) 
tempD[nf, 4] <- as.numeric(substr(tf[3], 48, 52)) 
# get right range of years 
year <- as.numeric(substr(tf, 1, 4)) 
use <- !is.na(year) & year < 2009 & year >1858 
year<- year[use] 
# extract year data and divide by 10 to convert to C 
tempD[nf,5:154] <- as.numeric(substr(tf[use], 116, 118))/10 
} 
Read 3560 items
[1] "homogenised" "mean"        "temperature" "2009-07-06"  "09:15"       "ABBR."      
Read 4698 items
Read 4014 items
Read 4885 items
Read 3957 items
Read 3292 items
Read 3460 items
Read 3276 items
Read 3293 items
Read 4737 items
# copy data to file so you dont have to extract again 
write.csv(tempD, file='tempD.csv',row.names=F) 
tempD <- read.csv(file='tempD.csv', header=T, as.is=T) 

  
 We pick off the name, longitude, latitude, and elevation in line 3. The last value in the line is the yearly temperature in 10ths of a degree Centigrade. We convert the temperatures to degree Centigrade before saving. The early years' temperatures, are missing(999999), as are the last months in 2009.
 6 The selected series 
  
print(tempD[,1:7]) 

 The first line is for Graz, longitude 15.45, latitude 47.08, altitude 377 metres, and then 9.0C, 7.7C,… for temperatures in 1859, 1860..
 This scraping was easy because the format was a 'flat' format. Each line of the data sets, except the first four, carried the same information about a different year. So we just needed to pick out the years we wished to include in the study. We have eaten our meat and vegetables, now for wine and dessert!
 7 Location of Stations 
 tiff('p/Station Locations.tif', w=880, h=680) 
boundary <- read.table(paste(dd, 'boundary.txt', sep='/'), header=F) 
tempD$lat[9] <- 46.55 
plot(boundary[,2], boundary[, 3], 
type='l',lwd=3, col='red', axes=F, xlab='', ylab='') 
points(tempD$long, tempD$lat, 
cex = tempD$alt/100, col='blue', pch=2, lwd=2, 
main="Long record Austrian temperature stations",cex=3) 
text(tempD$long, tempD$lat, tempD$loc, cex=1.6) 
dev.off(); picture('p/Station Locations.tif'  ) 

               
  The heights of stations are indicated by the blue triangles. Two of the stations are in the Alps.
 8 All stations yearly temperatures 
 year <- 1859:2008 
mt <- apply(tempD[, -(1:4)], 2, mean) 

 tiff('p/Annual Temperatures.tif', w=880, h=580) 
plot(year, tempD[1,-(1:4)], ylab='Yearly Average Temperature C', 
xlim=c(1800, 2008), ylim=c(-2, 11), cex=2 ) 
for( i in 1:10) { 
points(year, tempD[i, -(1:4)], col=i, pch=16) 
lines(year, tempD[i, -(1:4)], col=i) 
text(1800 + 5*i, tempD[i,5], substr(tempD$loc[i], 1, 6), 
cex=1.5, col=i) 
} 
title('Austrian Yearly Temperatures , 1859-2008', cex=2) 
text( 1850, 5, 'mean', cex=2) 
lines(year, mt, lwd=4) 
dev.off(); picture('p/Annual Temperatures.tif'  ) 

  We see a fairly similar pattern across all series, with the two mountain stations distinctly lower than the valley stations.
 9 Differences from overall mean 
 Subtract the mean for each year from the different station values for that year:
 rt <- t( t( tempD[,-(1:4)]) - mt) 

 tiff('p/Differences from mean.tif', w=880, h=580) 
plot(year, rt[1,] ,xlim=c(1800,2008), ylim=c(-8,6), 
ylab = 'Difference from overall mean', cex=2 ) 
for( i in 1:10) { 
points(year, rt[i, ], col=i, pch=16) 
lines(year, rt[i, ], col=i) 
text(1800+5*i, rt[i, 1], 
substr(tempD$loc[i], 1, 6), cex=1.5, col=i) 
} 
title('Austrian Temperatures Residuals, 1859-2008', cex=2) 
text( 1850, 0, ' mean', cex=2) 
lines(year, rep(0,length(year)), lwd=4) 
dev.off(); picture('p/Differences from mean.tif'  ) 

  Subtracting the mean from all stations much reduces the fluctuations, though the general shift difference between stations remains. Thus the temperature series for the different stations are at different levels, but change about the same each year.
 10  Adjust for yearly average and station average 
 Fit the additive model in which a temperature for a particular year and station is represented as a sum of terms consisting of a year value, a station value, and a residual:
 tiff('p/Residuals adjusting for year and stations.tif', w=880, h=580) 
ms <- apply(rt, 1, mean) 
rt <- rt - ms 
plot(year, rt[1,] ,xlim=c(1800, 2008), ylim=c(-1, 1), 
ylab= ' adjusted for station mean and yearly mean'  ) 
for( i in 1:10) { 
points(year, rt[i, ], col=i,pch=16) 
lines(year, rt[i, ], col=i) 
text(1800+5*i, rt[i, 1], substr(tempD$loc[i], 1, 6), 
cex=1.5, col=i) 
} 
title('Austrian Temperatures Res with station offsets, 1859-2008', cex=2) 
dev.off(); picture('p/Residuals adjusting for year and stations.tif'  ) 

  It looks like a mess, like a good residual plot should, nothing left to explain. But note that the curves are only varying by about 0.5C from the average curve, after adjusting for station offsets. This suggests that all the stations are following the same path, and will get a good idea of the climatology by just looking at the mean curve for all stations.
 11 Fitting curves to mean temperatures 
 We keep track of model fit with a version of the akaike information criterion:
AIC2 = 2 number of parameters - log likelihood.
 This criterion accepts an additional parameter if the likelihood ratio test rejects the null hypothesis omitting the parameter at the 5% significance level. 
 AIC2 <- function(s) 
2*s$df[1]+length(s$res)*log(max(1-s$r.sq,.00001))/2 
# subtract mean from square term to get orthogonal effects 
sqyear <- (year-mean(year))^2 
linsq <- lm(mt ~ year + sqyear) 
print(round( summary(linsq)$coef[-1, -(3:4)], 5)) 
cat('AIC2:', round(AIC2(summary(linsq)),1), '') 

 tiff('p/Fitting curves.tif', w=880, h=380) 
plot(year, mt ,xlim=c(1840, 2010), ylab ='Temperature C', 
ylim=c(5, 9) ,pch=16, col='green', 
main=' Parabola fitted to mean temperatures', cex=2) 
lines(year, mt) 
lines(year, linsq$fit, col='blue', lwd=3) 
dev.off(); picture('p/Fitting curves.tif'  ) 

  The zero AIC value occurs in fitting a single mean to all data points. The negative AIC here shows that the parabola, with two parameters, is superior to the single mean model.
 12 Parabola with lag 
 The lag model uses the previous year's temperature as one of the predicting variables in the model.
 mtlag1 <- c(mt[1], mt[ - length(mt)]) 
linsqlag <- lm(mt ~ year + sqyear + mtlag1) 
print(round( summary(linsqlag)$coef[-1, -(3:4)], 5)) 
cat('AIC2:', round(AIC2(summary(linsqlag)),1), '') 

  
 tiff('p/lag added.tif', w=880, h=480) 
plot(year, mt ,xlim=c(1840, 2010), 
ylim=c(5, 9) ,pch=16, col='green') 
lines(year, mt) 
lines(year, linsq$fit, col='blue', lwd=3) 
lines(year, linsqlag$fit, col='red', lwd=3) 
dev.off(); picture('p/lag added.tif'  ) 

  The additional lag term, which causes the wiggley fit, is not significant in the regression, and correspondingly produces an increase in the AIC2. Conclusion: drop the lag term.
 13 Fit a change point 
 For each year as a possible change point fit a different slope before and after the year, and then select that year with the smallest AIC2.
 options(digits=10) 
aic <- rep(0,146) 
for( change in 1861:2006){ 
# Construct a slope variable for the each segment, constant on other segments. 
y2 <- pmax(year - change, 0) 
y1 <- year - y2 
aic[change-1860]<-AIC2(summary( lm(mt ~ y1  + y2) )) 
} 
changept <- which( aic == min(aic) ) + 1860 
cat('AIC2:', round(min(aic), 1), '') 

  
 tiff('p/aic for different change points.tif', w=880, h=380) 
plot(1861:2006, aic, type='b', 
xlab='changepoint', ylab= 'AIC2', cex=2) 
title(' Change point for linear fit',cex=3) 
lines(1861:2006, rep(min(aic) + 2, length(aic)), col='red') 
text(1920, min(aic)+ 1, ' 95% confidence region', cex=2) 
arrows(1980, min(aic) + 1, 1980, min(aic)) 
dev.off(); picture('p/aic for different change points.tif'  ) 

  The red line in the plot corresponds to a decrease of 2 in the log likelihood; years where the plot lies below the redline give  a 95% confidence region for the change point. We see two possible minima in the curve, but the 95% interval is 1964 to 1989.  That's about where people say the sharp increase in temperature began.
 14 Fitted change point lines 
 y2 <- pmax(year - changept, 0) 
y1 <- year - y2 
reg <-  lm( mt ~ y1 +  y2 ) 
print(round(summary(reg)$coef[-1, -(3:4)], 5) ) 
cat('AIC2:', round(AIC2(summary(reg)), 1), '') 

 tiff('p/Fitted Change Point Lines.tif', w=880, h=480) 
plot(year, mt ,xlim=c(1840,2010), ylim=c(5, 9), 
ylab='Temperatures', xlab='Years', pch=16, col='green') 
lines( year, col = 'red', lwd = 2, reg$fit) 
lines( year, col = 'black', lwd = 1, mt) 
title(' Fitted lines with change point at 1980', cex=2) 
dev.off(); picture('p/Fitted Change Point Lines.tif'  ) 

  Not a bad fit, suggesting a rate of increase of .05 Centigrade per year since 1980. Note the improvement in the AIC2 compared to the parabola.
 15 Change in the individual stations 
 Find the best change point, at each station:
 # set up data frame to take change points and slopes 
chp <- data.frame( matrix(0, nrow=10, ncol=5)) 
names(chp) <- c('changept', 'slope1', 'se1', 'slope2', 'se2') 
rownames(chp) <- tempD$loc 
# run over different stations 
for (st in 1:10){ 
mt <- unlist(tempD[st, -(1:4)]) 
aic <- rep(0,146) 
for ( change in 1861:2006){ 
y2 <- pmax(year - change, 0) 
y1 <- year - y2 
aic[change-1860] <- summary(lm(mt ~ y1  + y2 ))$sigma 
} 
changept <- which(aic == min(aic)) + 1860 
y2 <- pmax(year - changept, 0) 
y1 <- year - y2 
 
reg <-  lm( mt ~ y1 +  y2 ) 
chp[st, 1] <- changept 
chp[st, 2:5] <- unlist(t(summary(reg)$coef[-1, -(3:4)])) 
} 
print(round(chp,4)) 

Good agreement, in change point estimates, and in slope after change, except for Vienna!
 16  What happened to Vienna 
 tiff('p/What happened to Vienna.tif', w=880, h=480) 
plot(1861:2006, aic, type='b', col='blue', ylab='AIC2') 
title('Vienna AIC2 for different change points', cex=2) 
dev.off(); picture('p/What happened to Vienna.tif'  ) 

  The optimal single breakpoint at 1886, was far earlier than for the other stations. Perhaps Vienna was growing faster then. The next local minimum occurs at 1968.
 17 Three segments for Vienna 
 Construct a slope variable for the each segment, constant on the other segments.
 y2 <- pmax(year - 1886, 0) 
y1 <- year - y2 
y3 <- pmax(year - 1968, 0) 
y23 <- y2 - y3 

 Regress with 2 segments. 
 reg2 <- lm(mt ~ y1  + y2 ) 
print(round(summary(reg2)$coef[-1, -(3:4)], 4) ) 
cat('AIC2:', round(AIC2(summary(reg2)), 1), '') 

 Regress with 3 segments
 reg3 <- lm(mt ~ y1 + y23 + y3) 
print(round(summary(reg3)$coef[-1, -(3:4)], 4) ) 
cat('AIC2:', round(AIC2(summary(reg3)), 1), '') 

 18 Three segment fit 
 tiff('p/ Two and Three Segments fitted.tif', w=880, h=680) 
plot(year, tempD[10,-(1:4)],pch=16,col='green') 
lines(year, tempD[10,-(1:4)]) 
lines(year,reg2$fit, lwd=3) 
lines(year,reg3$fit, lwd=3) 
title(' Fitting two and three segments to WIEN', cex=2) 
dev.off(); picture('p/ Two and Three Segments fitted.tif'  ) 

  A huge improvement in AIC2 for three segments vs 2; the suggested final slope is .35 which corresponds to .035 degrees centigrade per year. The three segment fit captures the rapidly declining values before 1880 without much affecting the estimate of slope of the recent rapid increase.
 19 Three segments for mean Austrian temperatures 
 Compute the AIC2 for for each choice of two change points , store in the matrix aic.
 mt <- apply(tempD[, -(1:4)], 2, mean) 
aic <- data.frame(matrix(0, 147, 147)) 
names(aic) <- 1860:2006 
row.names(aic) <- names(aic) 
for ( c1 in 1860:2005){ 
for (c2 in (c1+1):2006){ 
y2 <- pmax(year - c1, 0) 
y1 <- year - y2 
y3 <- pmax(year - c2, 0) 
y23 <- y2 - y3 
 
aic[c1-1859, c2-1859] <- AIC2(summary(lm(mt ~ y1 + y23 + y3))) 
} 
} 

 Look through the matrix for the pair of change points with the least aic.
 print(round(aic[17:23, 118:124],3)) 

  
 The change points are at 1879 and 1980 with the minimum AIC2 equal to -36.854.
 c1 <- 20+1859 
c2 <- 121+1859 
y2 <- pmax(year - c1, 0) 
y1 <- year - y2 
y3 <- pmax(year - c2, 0) 
y23 <- y2 - y3 
reg <- lm(mt ~ y1 + y23 + y3) 
print(round(summary(reg3)$coef[-1, -(3:4)], 4) ) 
cat('AIC2:', round(AIC2(summary(reg)), 1), '') 

 tiff('p/Fit with three segments.tif', w=880, h=480) 
plot(year, mt ,xlim=c(1840,2010),pch=16, col='green') 
lines(year, mt) 
lines( year, col = 'black', lwd = 3, lm(mt ~ y1 + y23 + y3)$fit  ) 
title(' Fitted lines with change points at 1879, 1980', cex=2) 
dev.off(); picture('p/Fit with three segments.tif'  ) 

  The three segment solution is a considerable improvement over the two segment solution. The two change point solution for the average curve  suggests an increase of 3.5C per century if the slope of the last segment continues (it won't) . The general story is of  decreasing temperatures up to 1879, a small rate of increase up to 1980, a large rate of increase since then.
 20 Algorithm to extrapolate using integrated brownian motion 
 Extrapolate the time series x using the integrated brownian motion model. This is a particular case of the Kalman filter. Assume x|z ~ N(z, v1) and diff(diff(z)) ~ N(0, v2). Here z is the fitted value to x, and is expected to be smoother than x. The second differences of z being independent normal characterises an integrated brownian motion. 

A straight line x would be fitted exactly by z under this model, with the v2 chosen to be zero. If x consisted of continuous straight line segments z would also closely approximated x. 

All the extrapolated values lie on the straight line through the last two fitted values.
 extrapolateIBM <- function(x, extra =0, iter=200){ 
# Extrapolate the time series x using integrated brownian motion model. 
# predictors for IBM ; begin with assumed small variance for IBM 
lx <- length(x) 
mx <- matrix(0, nrow=2*lx-2, ncol = lx) 
v1 <- mean(lm(x~seq(1:lx))$res^2) 
v2 <- v1/10 
 
# predictor matrix 
mx[1:lx, 1:lx] <- diag( rep(.0001, lx) ) 
for ( i in 1:(lx-2) )  mx[i+lx, i + 0:2] <- c(1, -2, 1) 
# Rotate x and z so that the rotated variables are independent. 
# svd holds for all v1 and v2 
s <- svd(mx) 
s$d[(lx-1):lx] <- 0 
d <- s$d 
vx <- t(s$v) %*% x 
# iterate using EM to estimate variances v1 and v2 
# the new estimate of v1 is the expected value of (vx-vz)^2, 
# the new estimate of v2 is the expected value of d^2 vz^2 , 
# when vz the posterior from vx~N(vz,v1), 0~N(vz*d , v2) 
# vx = norm( vx/(1 + d^2*v1/v2, 1/v1 + d^2/v2) 
d2 <- 1/d[-c(lx-1, lx)]^2 
vx2 <- vx[-c(lx-1, lx)]^2 
oldlike <- 0 
for(it in 1:iter){ 
vz <- diag( 1/(1 + d^2 * v1/v2) )  %*% vx 
vv <- 1/(1/v1 + d^2/v2) 
mv <- v1 + v2*d2 
loglike <- -0.5*( sum(log(mv) + vx2/mv) ) 
if(it==1 |loglike > oldlike +.001){ 
oldlike <- loglike 
#EM updates; accelerate a little 
acc <- 0 
if(it>5) acc <- 0.75 
if(it>10) acc <- 1 
v2 <-((1+acc)*sqrt(sum( d^2*vx^2*(vv/v1)^2+d^2*vv )/(lx-2))- 
acc*sqrt(v2) )^2 
vv <- 1/(1/v1 + d^2/v2) 
vz <- diag( 1/(1 + d^2 * v1/v2) )  %*% vx 
 
v1 <- ( (1+acc)*sqrt(mean( vx^2 * ( 1- vv/v1 )^2  + vv )) - acc * sqrt(v1) )^2 
 
} 
} 
# two additional parameters over fitting a mean and a variance 
# the third term is minus log like for fitting a mean and variance 
aic2 <- 4 - oldlike + length(x) * (log(var(x)) + 1)/2 
z <- s$v %*% diag( 1/(1 + d^2 * v1/v2) )  %*% vz 
# do the regression with the given v1, v2 to get estimates 
# and compute se for the present and future z's 
mx[1:lx, 1:lx] <- diag(1, lx) 
w <- rep(1/v1, 2*(lx-1)) 
w[(lx+1) : (2*(lx-1))] <- 1/v2 
 
xx <- c(x, rep(0, lx-2)) 
cov <- summary(lm(xx ~ -1+mx, weights=w))$cov 
se <- sqrt(diag(cov)) 
# extrapolate with se 
k <- 1:extra 
zz <- z[lx] + k * (z[lx] - z[lx-1]) 
see <- sqrt((k+1)^2*cov[lx,lx] - 2*k*(k+1)*cov[lx, lx-1] + (k+1)^2*cov[lx-1, lx-1]) 
 
return(list( est=c(z,zz), se=c(se,see), aic2 = aic2 )) 
} 

 21 Extrapolate the next twenty years, until 2028 
 temp <- read.csv('tempD.csv', header=T, as.is=T) 
data <- tempD[, -(1:4)] 
mt <- apply(data, 2, mean) 
f<-extrapolateIBM(mt, extra=20) 
print(round(f$aic2, 1)) 

 tiff('p/Extrapolate next twenty years.tif', w=880, h=580) 
plot(1859:2028, f$est, type='l', ylim= c(5,12), ylab='Austrian Temp in Degrees C', xlab='years', lwd=3 ) 
lines(1859:2028, f$est) 
lines(1859:2028, f$est + 2*f$se, col='red', lwd=2) 
lines(1859:2028, f$est - 2*f$se, col='blue', lwd=2) 
for(y in 5:12) lines(1859:2028, rep(y,170), col='grey') 
title(' Integrated brownian motion extrapolation +- 2se', cex=3) 
points(1859:2008, mt, pch=16) 
dev.off(); picture('p/Extrapolate next twenty years.tif'  ) 

  I like this better than the 3 segment solution, which assesses the break points too confidently, and suggests too strongly that the slope of the last segment is the one to be sustained. After a small initial decrease, the smoothed curve increases over the last 100 years, and more rapidly in the last 30 years. The suggested future increase is 4C in the next 100 years, with a very wide error band about that prediction. Similar extrapolations occur for the individual stations.
 22 Conclusion 
 Austria is warming up. The worldwide data suggests that the 30-50 degrees latitude region, full of population, production, and pollution, are producing lots of carbond dioxide, and warming up faster than other latitudes.