rawData <- "Year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec J-D D-N DJF MAM JJA SON Year 1881 -26 -30 -6 -10 -10 -35 -18 -15 -29 -33 -42 -34 -24 -24 -28 -9 -23 -35 1881 1882 -16 -5 -13 -33 -30 -38 -32 -19 -28 -42 -42 -57 -30 -28 -18 -26 -29 -37 1882 1883 -51 -45 -22 -27 -32 -15 -14 -25 -35 -37 -37 -33 -31 -33 -51 -27 -18 -37 1883 1884 -31 -23 -39 -42 -40 -40 -32 -22 -35 -37 -39 -34 -35 -34 -29 -40 -31 -37 1884 1885 -65 -36 -27 -44 -37 -45 -33 -30 -26 -27 -23 -9 -34 -36 -45 -36 -36 -25 1885 1886 -39 -46 -31 -24 -22 -35 -11 -25 -23 -34 -36 -32 -30 -28 -32 -26 -24 -31 1886 1887 -66 -53 -37 -44 -31 -26 -11 -29 -27 -40 -34 -42 -37 -36 -50 -37 -22 -34 1887 1888 -50 -49 -47 -37 -30 -26 -22 -25 -21 -13 -8 -24 -29 -31 -47 -38 -24 -14 1888 1889 -25 6 -4 -5 -9 -15 -19 -27 -23 -33 -40 -36 -19 -18 -14 -6 -20 -32 1889 1890 -54 -43 -37 -39 -52 -41 -33 -40 -44 -26 -53 -37 -42 -41 -45 -43 -38 -41 1890 1891 -52 -57 -20 -33 -23 -26 -26 -21 -19 -26 -41 -12 -30 -32 -48 -26 -24 -29 1891 1892 -39 -13 -37 -45 -34 -24 -33 -30 -22 -31 -54 -50 -34 -31 -21 -38 -29 -36 1892 1893 -88 -61 -21 -39 -41 -31 -11 -27 -25 -19 -20 -35 -35 -36 -67 -34 -23 -21 1893 1894 -53 -34 -27 -46 -40 -46 -23 -26 -36 -30 -41 -31 -36 -36 -41 -38 -32 -35 1894 1895 -57 -51 -33 -28 -31 -23 -22 -20 -13 -18 -15 -22 -28 -29 -47 -31 -22 -15 1895 1896 -28 -22 -34 -40 -19 -16 -10 -16 -13 -4 -18 -15 -20 -20 -24 -31 -14 -11 1896 1897 -25 -22 -22 -10 -7 -15 -8 -11 -17 -15 -23 -20 -16 -16 -21 -13 -11 -18 1897 1898 -10 -30 -52 -29 -39 -22 -24 -25 -25 -33 -43 -31 -30 -29 -20 -40 -24 -34 1898 1899 -26 -40 -33 -19 -19 -31 -13 -10 -11 -5 8 -32 -19 -19 -33 -24 -18 -3 1899 1900 -38 2 -4 -13 -10 -12 -13 -13 -8 -3 -17 -11 -12 -13 -23 -9 -13 -9 1900 1901 -26 -5 -2 -7 -19 -15 -20 -21 -23 -29 -26 -29 -18 -17 -14 -9 -19 -26 1901 1902 -17 2 -26 -28 -32 -33 -21 -30 -29 -35 -43 -47 -28 -27 -14 -29 -28 -36 1902 1903 -27 4 -17 -37 -34 -44 -29 -39 -45 -43 -34 -46 -33 -33 -24 -29 -37 -40 1903 1904 -57 -47 -37 -45 -43 -34 -35 -35 -41 -33 -14 -22 -37 -39 -50 -42 -35 -29 1904 1905 -31 -56 -23 -36 -31 -29 -25 -20 -20 -28 -12 -23 -28 -28 -36 -30 -25 -20 1905 1906 -29 -33 -19 -6 -22 -15 -24 -14 -24 -18 -41 -17 -22 -22 -28 -16 -18 -28 1906 1907 -46 -52 -30 -49 -52 -47 -41 -42 -31 -27 -45 -46 -42 -40 -38 -44 -44 -35 1907 1908 -44 -29 -49 -44 -37 -29 -27 -37 -25 -33 -42 -41 -36 -37 -39 -43 -31 -33 1908 1909 -58 -40 -45 -47 -44 -43 -35 -22 -24 -24 -20 -41 -37 -37 -47 -45 -33 -23 1909 1910 -36 -35 -38 -37 -38 -37 -24 -26 -32 -35 -44 -59 -37 -35 -37 -38 -29 -37 1910 1911 -55 -51 -54 -52 -45 -43 -28 -31 -28 -22 -21 -22 -38 -41 -55 -50 -34 -23 1911 1912 -26 -13 -34 -25 -23 -23 -42 -55 -47 -55 -35 -38 -35 -33 -20 -27 -40 -46 1912 1913 -42 -37 -41 -40 -42 -48 -33 -29 -30 -35 -23 -6 -34 -36 -39 -41 -36 -29 1913 1914 -6 -13 -23 -27 -19 -21 -21 -18 -20 -3 -21 -18 -18 -17 -8 -23 -20 -15 1914 1915 -24 -9 -16 -2 -5 -4 -3 -13 -8 -22 -13 -21 -12 -11 -17 -8 -7 -14 1915 1916 -19 -18 -34 -28 -34 -39 -29 -22 -28 -21 -38 -64 -31 -28 -19 -32 -30 -29 1916 1917 -42 -51 -52 -40 -57 -38 -19 -22 -17 -38 -31 -71 -40 -39 -52 -50 -26 -29 1917 1918 -51 -45 -30 -47 -48 -37 -30 -40 -26 -10 -23 -36 -35 -38 -56 -42 -36 -20 1918 1919 -19 -23 -21 -13 -20 -21 -22 -17 -14 -20 -46 -35 -23 -23 -26 -18 -20 -27 1919 1920 -14 -23 -2 -16 -14 -22 -28 -22 -23 -29 -32 -47 -23 -22 -24 -11 -24 -28 1920 1921 -5 -18 -17 -20 -23 -17 -3 -29 -22 -11 -19 -20 -17 -19 -23 -20 -16 -18 1921 1922 -33 -39 -22 -28 -33 -28 -18 -29 -32 -29 -19 -21 -28 -27 -31 -28 -25 -27 1922 1923 -27 -29 -26 -38 -39 -24 -30 -32 -29 -11 0 -5 -24 -25 -26 -34 -29 -13 1923 1924 -23 -27 -13 -35 -24 -23 -21 -24 -27 -27 -16 -34 -25 -22 -19 -24 -23 -23 1924 1925 -33 -28 -20 -21 -27 -31 -24 -17 -20 -24 -4 12 -20 -24 -32 -22 -24 -16 1925 1926 17 10 20 -12 -18 -19 -13 -4 -8 -3 -9 -22 -5 -2 13 -3 -12 -7 1926 1927 -21 -8 -32 -27 -24 -21 -11 -19 -11 0 -4 -33 -17 -17 -17 -28 -17 -5 1927 1928 -1 -5 -22 -28 -28 -31 -14 -19 -14 -15 -3 -11 -16 -18 -13 -26 -21 -11 1928 1929 -35 -46 -27 -34 -34 -36 -30 -27 -23 -11 -6 -46 -30 -27 -31 -32 -31 -13 1929 1930 -20 -18 -6 -20 -22 -18 -13 -7 -12 -10 13 -6 -12 -15 -28 -16 -13 -3 1930 1931 -4 -20 -4 -16 -13 0 8 0 -10 2 -5 -2 -5 -6 -10 -11 3 -4 1931 1932 13 -16 -15 -6 -15 -23 -15 -13 -1 -1 -19 -19 -11 -10 -2 -12 -17 -7 1932 1933 -26 -23 -23 -20 -20 -23 -11 -14 -22 -18 -30 -44 -23 -21 -23 -21 -16 -23 1933 1934 -20 -3 -29 -28 -10 -11 -6 -8 -18 -5 5 -1 -11 -15 -22 -22 -8 -6 1934 1935 -31 19 -8 -30 -28 -19 -15 -16 -17 -3 -26 -16 -16 -15 -4 -22 -17 -15 1935 1936 -23 -31 -17 -14 -13 -10 7 -3 -5 1 4 7 -8 -10 -24 -15 -2 0 1936 1937 -1 16 -10 -12 -1 1 5 6 16 12 15 -5 4 5 7 -8 4 14 1937 1938 12 6 20 18 5 -5 4 3 10 17 13 -12 8 8 4 14 1 13 1938 1939 -3 -5 -18 -8 -5 -3 1 -6 -3 -7 2 34 -2 -6 -7 -10 -3 -3 1939 1940 -20 -2 2 10 2 1 15 -5 5 0 5 9 2 4 4 5 4 3 1940 1941 7 19 1 12 7 4 9 -1 -15 17 9 18 7 7 12 6 4 4 1941 1942 19 -10 -3 -2 5 2 -1 -3 -1 0 0 2 1 2 9 0 -1 0 1942 1943 -15 9 -9 6 14 0 13 4 2 21 16 22 7 5 -1 4 6 13 1943 1944 33 22 19 9 18 10 17 14 24 22 14 10 18 19 26 15 14 20 1944 1945 14 10 11 16 0 -3 1 12 -2 6 4 -11 5 7 11 9 3 3 1945 1946 11 5 -5 4 -9 -17 -5 -21 -5 -6 -6 -29 -7 -6 2 -3 -14 -6 1946 1947 -6 -3 13 4 -4 -9 -5 -6 -8 11 10 -13 -1 -3 -13 4 -7 4 1947 1948 11 -11 -10 -6 3 3 -11 -9 -9 -1 -6 -22 -6 -5 -4 -4 -6 -5 1948 1949 13 -18 -2 -10 -10 -22 -11 -6 -7 -2 -7 -11 -8 -9 -9 -7 -13 -5 1949 1950 -26 -29 -3 -19 -15 -7 -10 -20 -13 -18 -30 -12 -17 -17 -22 -12 -13 -20 1950 1951 -33 -41 -19 -10 -1 -6 -3 9 7 12 1 15 -6 -8 -29 -10 0 7 1951 1952 13 12 -9 4 -2 1 5 8 6 -6 -16 -5 1 3 13 -3 5 -5 1952 1953 9 16 14 17 6 2 3 8 8 8 0 13 9 7 7 12 4 5 1953 1954 -19 -5 -10 -11 -18 -14 -23 -19 -14 -4 11 -13 -12 -10 -4 -13 -19 -3 1954 1955 17 -15 -31 -18 -20 -14 -8 10 -9 0 -24 -29 -12 -10 -4 -23 -4 -11 1955 1956 -17 -26 -23 -22 -25 -16 -10 -27 -16 -22 -18 -9 -19 -21 -24 -24 -18 -19 1956 1957 -7 2 1 9 10 18 3 15 7 4 13 17 8 6 -4 7 12 8 1957 1958 42 23 15 4 7 -9 7 -2 -4 4 5 4 8 9 27 8 -1 2 1958 1959 11 10 21 17 9 6 6 0 -8 -4 -8 3 5 5 9 16 4 -7 1959 1960 0 18 -31 -14 -6 -4 -2 5 5 5 -9 19 -1 -3 7 -17 0 0 1960 1961 4 20 10 11 21 13 0 3 7 6 5 -15 7 10 14 14 5 6 1961 1962 5 17 12 11 -10 8 -4 -6 3 0 7 0 4 2 2 4 0 3 1962 1963 1 20 -13 -8 -3 6 14 26 23 9 14 3 8 7 7 -8 15 15 1963 1964 -4 -9 -27 -33 -27 -4 -5 -24 -37 -29 -20 -30 -21 -18 -3 -29 -11 -29 1964 1965 -10 -18 -9 -19 -6 -11 -21 -6 -18 -7 -8 -6 -12 -14 -19 -12 -13 -11 1965 1966 -17 -1 10 -11 -6 -2 10 -5 1 -15 0 -3 -3 -3 -8 -3 1 -4 1966 1967 -6 -24 8 -2 12 -6 6 1 0 11 -1 -1 0 0 -11 6 0 3 1967 1968 -22 -14 26 -4 -9 -1 -7 -4 -13 13 -2 -12 -4 -3 -12 4 -4 -1 1968 1969 -6 -8 -2 19 13 13 -3 0 9 15 16 31 8 4 -9 10 3 13 1969 1970 9 23 10 4 -4 -1 -1 -13 12 4 4 -12 3 6 21 3 -5 7 1970 1971 -2 -20 -21 -11 -10 -20 -11 -3 -1 -4 -7 -11 -10 -10 -11 -14 -11 -4 1971 1972 -25 -20 -3 -1 -2 6 0 21 5 5 -1 19 0 -2 -19 -2 9 3 1972 1973 25 29 26 25 24 17 11 2 8 13 4 -8 15 17 24 25 10 8 1973 1974 -14 -25 -4 -11 -3 -4 -1 12 -8 -6 -10 -10 -7 -7 -15 -6 2 -8 1974 1975 3 2 12 -1 20 -1 -2 -20 -5 -9 -15 -20 -3 -2 -2 10 -8 -9 1975 1976 -6 -9 -26 -14 -28 -14 -12 -18 -9 -27 -10 1 -14 -16 -12 -23 -14 -16 1976 1977 11 16 17 21 30 24 21 18 -2 -3 15 4 14 14 10 23 21 3 1977 1978 6 10 14 10 3 -6 7 -19 8 -1 9 2 4 4 7 9 -6 5 1978 1979 8 -16 11 12 -3 6 -5 11 21 20 21 41 11 7 -2 7 4 21 1979 1980 23 31 24 27 29 11 23 19 15 11 23 11 21 23 32 27 18 16 1980 1981 50 38 45 25 17 23 30 33 14 7 19 30 28 26 33 29 29 13 1981 1982 2 10 -9 0 13 2 14 -2 3 4 4 37 7 6 14 1 5 4 1982 1983 47 36 38 28 31 17 12 31 37 11 23 11 27 29 40 33 20 23 1983 1984 24 10 25 5 31 0 13 12 18 5 -4 -13 10 12 15 20 8 6 1984 1985 18 -10 12 5 7 15 -4 11 7 7 2 10 7 5 -2 8 7 5 1985 1986 23 38 26 21 17 7 7 8 0 7 3 8 14 14 24 21 7 3 1986 1987 27 39 12 18 17 33 42 20 31 27 22 45 28 25 25 16 32 27 1987 1988 52 37 47 36 39 40 28 33 31 30 0 21 33 35 45 41 33 20 1988 1989 7 30 30 20 8 8 30 27 33 28 10 30 22 21 19 19 22 23 1989 1990 31 30 69 47 38 29 44 25 17 39 38 36 37 36 30 51 33 31 1990 1991 35 46 30 45 32 50 49 40 41 21 20 26 36 37 39 36 46 27 1991 1992 39 35 40 17 24 18 2 1 -10 -4 -6 11 14 15 33 27 7 -7 1992 1993 29 30 29 17 18 10 14 3 0 15 0 13 15 15 23 21 9 5 1993 1994 23 -6 20 31 17 36 23 16 32 37 40 29 25 24 10 23 25 37 1994 1995 44 72 42 42 21 38 50 40 24 44 38 25 40 40 49 35 42 35 1995 1996 25 48 32 27 21 19 36 49 25 17 35 35 31 30 33 26 35 26 1996 1997 28 33 49 33 32 53 29 36 44 53 58 58 42 40 32 38 39 52 1997 1998 55 84 60 57 68 73 66 65 44 41 45 52 59 60 65 62 68 43 1998 1999 42 60 28 26 23 36 31 30 31 34 30 35 34 35 51 26 32 32 1999 2000 18 53 51 54 31 38 37 40 36 20 29 25 36 37 35 45 38 28 2000 2001 39 41 55 44 53 46 54 48 51 47 66 51 50 47 35 51 49 55 2001 2002 71 71 89 56 57 49 58 47 54 51 51 38 58 59 64 67 51 52 2002 2003 68 51 54 49 53 42 50 65 61 68 50 69 57 54 52 52 52 60 2003 2004 54 67 59 52 36 35 21 43 48 60 67 47 49 51 63 49 33 58 2004 2005 69 55 69 63 55 59 55 58 68 73 64 64 63 61 57 62 57 68 2005 2006 48 63 59 45 42 55 43 64 55 60 65 72 56 55 58 48 54 60 2006 2007 89 64 65 68 62 54 56 57 52 55 49 40 59 62 75 65 56 52 2007 2008 16 26 65 44 41 35 54 36 53 56 58 49 44 44 28 50 42 56 2008 2009 55 45 48 49 54 62 67 56 66 60 66 60 57 56 50 50 62 64 2009 2010 69 75 85 77 65 56 51 55 54 63 72 45 64 65 68 76 54 63 2010 2011 46 44 57 56 43 51 66 66 50 55 47 43 52 52 45 52 61 50 2011" library(ggplot2) library(car) library(lmtest) set.seed(1) #setwd("~/Research/blog/climate") theData <- read.table(textConnection(rawData),header=TRUE,na.string="*****") # There has to be a more elegant way to do this theData$means = rowMeans(aggregate(theData[,c("DJF","MAM","JJA","SON")], by=list(theData$Year), FUN="mean")[,2:5]) #4th order autoregressive is best ar(dat$means) rawChanges = diff(theData$means, 1) # the real data #rawChanges = diff(sim, 1) # the simulated data # SD on yearly changes sd(rawChanges) # Subtract off the mean, so that the distribution now has an expectaion of zero changes = rawChanges - mean(rawChanges) png("randomWalk.png") # First simulation, with plotting plot.ts(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3, xlab="Year", ylab="Temperature anomaly in hundreths of a degrees Celsius") trials = 1000 finalResults = rep(0,trials) for(i in 1:trials) { jumps = sample(changes, 130, replace=T) # Add lines to plot for this, note the "alpha" term for transparency lines(cumsum(c(0,jumps)), col=rgb(0, 0, 1, alpha = .1)) finalResults[i] = sum(jumps) } # Re-plot red line again on top, so it's visible again lines(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3) dev.off() png("meanReg.png") qplot(rawChanges[-length(rawChanges)],rawChanges[-1]) + geom_smooth(method="lm") + theme_bw() + xlab("YoY Change Last Year") + ylab("YoY Change This Year") + ggtitle("Classic regression to the mean") dev.off() dat <- data.frame(means=theData$means,year=theData$Year) #fit a simple curve to the data png("means.png") ggplot(dat,aes(x=year,y=means)) + geom_point() + geom_smooth(method="lm",formula=y~poly(x,3)) + theme_bw() dev.off() model <- lm(formula=means ~ poly(year,3),data=dat,na.action=na.omit) Anova(model,type='II') residSd <- sd(resid(model)) dat$resid <- resid(model) #is there a trend left in the residuals? png("resid.png") qplot(year,resid,data=dat)+theme_bw() dev.off() #Is there year to year correlation in the residuals? #acf(dat$resid) dwtest(model) png("simulated.png") sim <- predict(model) + rnorm(nrow(dat),sd=residSd) ggplot(dat,aes(x=year,y=sim)) + geom_point() + theme_bw() dev.off() model <- lm(formula=sim ~ poly(year,3),data=dat,na.action=na.omit) Anova(model,type='II') # # # now we are going to run the analysis from # http://www.statisticsblog.com/2012/12/the-surprisingly-weak-case-for-global-warming/ # using simulated data from our linear model (which was very signficant) # # # Get a single vector of Year over Year changes #rawChanges = diff(theData$means, 1) # the real data rawChanges = diff(sim, 1) # the simulated data # SD on yearly changes sd(rawChanges) # Subtract off the mean, so that the distribution now has an expectaion of zero changes = rawChanges - mean(rawChanges) # Find the total range, 1881 to 2011 (theData$means[131] - theData$means[1])/100 # Year 1 average, year 131 average, difference between them in hundreths y1a = theData$means[1]/100 + 14 y131a = theData$means[131]/100 + 14 netChange = (y131a - y1a)*100 png("randomWalkSim.png") # First simulation, with plotting plot.ts(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3, xlab="Year", ylab="Temperature anomaly in hundreths of a degrees Celsius") trials = 1000 finalResults = rep(0,trials) for(i in 1:trials) { jumps = sample(changes, 130, replace=T) # Add lines to plot for this, note the "alpha" term for transparency lines(cumsum(c(0,jumps)), col=rgb(0, 0, 1, alpha = .1)) finalResults[i] = sum(jumps) } # Re-plot red line again on top, so it's visible again lines(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3) dev.off() # Re-plot red line again on top, so it's visible again lines(cumsum(c(0,rawChanges)), col="red", ylim=c(-300,300), lwd=3) ( length(finalResults[finalResults>sum(rawChanges)]) + length(finalResults[finalResults<(-sum(rawChanges))]) ) / trials # Plot of YoY changes over time plot(rawChanges,pch=20,col="blue", xlab="Year", ylab="YoY change (in hundreths of a degree)") # Is there a trend? absRawChanges = abs(rawChanges) pts = 1:130 summary(lm(absRawChanges~pts)) png("meanRegSim.png") qplot(rawChanges[-length(rawChanges)],rawChanges[-1]) + geom_smooth(method="lm") + theme_bw() + xlab("YoY Change Last Year") + ylab("YoY Change This Year") + ggtitle("Classic regression to the mean") dev.off()