Op dezelfde manier als in de vorige post analyseer ik de temperaturen in Nederland. Het aardige is dat er een maandelijkse historie is van de gemiddelde dagtemperatuur in De Bilt, een serie die terug gaat tot 1706. Om een bruikbaar plaatje te maken, moeten we eerst de seizoensvariatie neutraliseren. We moeten dus eerst de best passende cyclische functie vinden, bijvoorbeeld geschikte $$a$$ en $$b$$ voor $$a*sin(2*pi*T+b)$$, waarbij $$T$$ het jaar is. Herinner je de formule voor $$sin(a+b)$$. We kunnen dus een lineaire regressie maken op de sinus en cosinus. Voor de zekerheid nemen we ook nog een dubbele periode mee, dus maken we een regressie met 4 parameters voor de periode: $$sin(2*pi*T)$$, $$sin(4*pi*T)$$, $$ cos(2*pi*T)$$ en $$cos(4*pi*T)$$. 
Halen we de seizoensinvloed weg, dan krijgen we een “genormaliseerde” reeks. Dat is het blauw op de achtergrond. Nog steeds zijn de variaties fors, wel 10 keer zo groot als de gemiddelden op de globale grafieken. Dat geeft verder niet echt. Verder heb ik gekeken of in dit signaal ook de zonnecyclus, de variaties van El Nino en de gevolgen van vulkanische uitbarstingen te zien zijn. Ze zijn er inderdaad alle drie, maar het verband met El Nino (ENSO) is toevallige ruis, die heb ik weggelaten. De effecten van vulkanen zijn wel duidelijk maar pas laat, na 20 maanden, goed zichtbaar. In het plaatje hiernaast zie je dat de statistiek maximaal 0.4 graad van de variatie in opwarming toerekent aan de zonnecyclus en dat vulkanen dit beeld de andere kant uit met ongeveer 0.2 graad tijdelijk kunnen verstoren. Het gehele beeld is een gemiddelde verhoging van de temperatuur met 0.39 per decennium, ofwel 1.7 graden Celcius sinds 1970.
Code
# haal gecachte data op vanuit mijn bibliotheek trial <- read.lab() # Labrije tijdsreeks voor De Bilt trial <- window(trial, start=1970) read.mei() # ENSO read.ssp() # sun spots read.vol() # Vulkanen pvol=20 pssp=1 pmei=5 vol.lag <- lag(vol.ts, -pvol) ssp.lag <- lag(ssp.ts, -pssp) mei.lag <- lag(mei.ts, -pmei) start.y <- 1970 end.y <- end(trial) debilt <- window(trial, start=start.y, end = end.y) mei <- window(mei.lag, start=start.y, end = end.y) ssp <- window(ssp.lag, start=start.y, end = end.y) vol <- window(vol.lag, start=start.y, end = end.y) year <- index(trial) f1 <- sin(2*pi*year) f2 <- cos(2*pi*year) f3 <- sin(4*pi*year) f4 <- cos(4*pi*year) fanta <- lmRob(debilt~year+f1+f2+f3+f4+ssp+vol+mei) # Robust regression periodic.fit <- f1*fanta$coefficients['f1'] + f2*fanta$coefficients['f2'] + f3*fanta$coefficients['f3'] + f4*fanta$coefficients['f4'] ssp.fit <- ssp*fanta$coefficients['ssp'] vol.fit <- vol*fanta$coefficients['vol'] mei.fit <- mei*fanta$coefficients['mei'] year.fit <- fanta$fit - (ssp.fit+periodic.fit+vol.fit+mei.fit)
Andere data en een vergelijkbare schatting zijn te vinden op het Berkely Earth project: http://berkeleyearth.lbl.gov/regions/netherlands