I saw that GISTemp had updated, and decided to take a breather from the sewing spree. For those wondering, here are the current uncertainty bounds on temperature trends since 2000.
Update 3:36 pm. Due to the issue with reading GIStemp automatically, I had to modify my script, I modified in one spot but not another. Consquently, I think the graphs below are really NOAA/NCDC. I’ll be updating tomorrow and adding in the models.
Because the SRES for the AR4 were not frozen in Jan 2000, I prefer to do comparisons since Jan 2001. Here is the GISTemp since that time:
Note: My “red” uncertainty intervals use closed form correction to Quenoulli which attempts to widen the 95% uncertainty intervals to account for the known excess false positives of that method. (Unfortunately, the notes in my script don’t mention the paper– but I could dig that up.) The ARIMA(4,0,0) is the ARIMA(p,0,q) with p≤4 and q≤4 that happens to give the widest uncertainty intervals of all possibilities that are stationary.
Currently, the linear trends since 2000 or 2001 are not statistically significant based if we assume ‘noise’ can be described using the noise models provided and we use the 2-sigma uncertainty intervals as our cut-off point for significance. Meanwhile, a nominal trend of 0.2C/decade continues to fall outside the 2 sigma range of possible trends and would be rejected as inconsistent with the observation using these methods.
Bleg
I’d show models and compare to other agencies but Hadcrut4 has not reported and — more importantly– I want to ask those who have run R if they have found any issues querying “http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt” using R. That page returned a 403 when I tried to use R. Of course, I immediately tried my browser, which worked fine. So I tried spoofing the user agent, and was rejected at the GISTemp page. Does anyone know if there are changes in the access policy? Could some magical incantation permit me to access? (BTW: As someone who seriously checks all access to my blog to inhibit bots, my question is not an attempt to complain about GISS, I’m just trying to discover if there has been some change in access policy. If I have to manually download and read a copy off my hard disk, that’s ok. But I want to know because having the script grab the data is more convenient and if others can still download using R, it may just be some problem with headers I’m passing or something.)
403 forbidden.
Lucia,
Yes, I noted the same problem here. I run a daily wget to check on updates to the main indices, and since the GISS server troubles, they haven’t responded. The others all do. My guess is that when the servers ran into trouble, someone disabled automated downloads etc to lighten the load. Hopefully temporarily.
Could this be connected to this messsage which is currently being displayed on the NASA/GISS website:
— Please Note —
Due to technical problems with the GISS webserver, most interactive content (creating scientific plots using web forms, searching the publications database, etc.) is currently disabled.
Ray–
It shouldn’t be related. In principle, if I can load the web page with my browser, I should be able to read it. I would image that particular page is static– just updated every now and then. That said, it’s not certain.
I wonder if it does lighten the load? If the page is static, it won’t make much of a difference. (On the other hand, I am a person who is blocking scrapers from racing through my images. Those are static, and I still don’t like things scraping mercilessly!)
At least I know someone else is having the same issue. I worked around it by loading manually and saving the text to a file manually. It’s not a big deal— but obviously, I’d rather just run the script and have it read the page!
Lucia
Their server seems to be working again but it is difficult to get the original ‘raw data’ archived station data. I found that individual stations I bookmarked do not work and so had to ‘reboookmark’ these.
There seems to be a difference between ‘after removing suspicious records’ and ‘after GISS homogeneity adjustment’.
eg De Bilt
V1 (original raw data – now archived)
http://data.giss.nasa.gov/cgi-bin/gistemp/show_station.cgi?id=633062600003&dt=1&ds=1
V3
http://data.giss.nasa.gov/cgi-bin/gistemp/show_station.cgi?id=633062600000&dt=1&ds=13
V4
http://data.giss.nasa.gov/cgi-bin/gistemp/show_station.cgi?id=633062600000&dt=1&ds=13
I would say that V1, the original ‘raw data’, will disappear soon.
Some of the US w/s (eg Boulder) seem to show the ‘rawer temps’ at V3 and their adjustment at V4.
GISS has 2012 at 0.55C above the 1951-1980 mean and Jan 2013 at 0.62C.
Jan 2013
http://data.giss.nasa.gov/cgi-bin/gistemp/nmaps.cgi?year_last=2013&month_last=1&sat=4&sst=3&type=anoms&mean_gen=01&year1=2013&year2=2013&base1=1951&base2=1980&radius=1200&pol=reg
Ian George,
“GISS has 2012 at 0.55C above the 1951-1980 mean and Jan 2013 at 0.62C.”
Any idea why that map shows 0.62c for January when the data files show 0.61c?
The annual figures are both 0.55c
No idea. But it could be rounding. It could also be that data are trickling in and values were computed at slightly different times based on slightly different data. You’d have to ask someone at GISS to get informed guesses.
Lucia,
“Currently, the linear trends since 2000 or 2001 are not statistically significant..”
I assume you mean “not significantly different from zero trend” at 95% confidence. Or do you mean, “not statistically different from each other” at 95% confidence?
Looks like GISS is blocking HTTP 1.0 requests. I’m unable to download with R or wget, but can with curl. Using curl in HTTP 1.0 mode fails, so it’s not a user agent check. It doesn’t appear to be a vhost issue, as the Host header is sent by R.
There’s a R package for curl (http://cran.r-project.org/web/packages/RCurl/index.html)
lucia. post your code for downloading the url.
there are a variety of ways to do it.
download.file() or you can get fancier with RCurl and get diagnostics on your requests.. and submit forms and passswords etc.
Yes. That’s what I mean.
I thought it might be the http 1.0 issue; I figured it can’t be a user agent block because then spoofing firefox would have fixed the problem. I wanted to find out if others are having a problem before I started trying to tease out the issue.
Mosher–
Agency=”GISTemp”
dataSource=”http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt”
dataSource=”http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt”
#dataSource=”http://data.giss.nasa.gov/gistemp/_tabledata3/GLB.Ts+dSST.txt”
# http://data.giss.nasa.gov/gistemp/_graphs3/Fig.C.txt
# http://data.giss.nasa.gov/gistemp/_graphs3/Fig.C.txt
url <- c("http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt") url <- c("http://data.giss.nasa.gov/gistemp/_tabledata3/GLB.Ts+dSST.txt") url <- c("http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt") file <- c("GLB.Ts+dSST.txt") download.file(url, file) As you can see... I leave in cruft....
mt (Comment #110553)
Thanks, mt. I switched from wget to cURL, and it worked. cURL also allows me to download only where there is a new version, so I think I’ll stay with them.
Can we have a 30 year trend?
bugs–
The data are at http://data.giss.nasa.gov/gistemp/tabledata_v3/GLB.Ts+dSST.txt . “We” can grab it and have any length trend “we” want. “We” includes you. Knock yourself out.
bugs (#110559): “Can we have a 30 year trend?”
Here are around 5000 of them. Enjoy.
HaroldW,
Nice graph.
SteveF/HaroldW
I suspect it’s not the graph bugs hoped to see.
You say red is the “corrected” OLS trend. What correction do you use?
I’m asking because I’m getting a regression slope of +0.00132 for GISS data during the same time period (Jan 2001 thru Jan 2013) using standard OLS.
lucia, trying adding mode = “wb”
download.file(url, mode=”wb”)
I’ll try in a little bit
Lucia,
No, probably not. But still a very nice graph; nice in the sense that it provides a wealth of information that is 100% consistent with significant pseudo-cyclical behavior over the instrument era.
Kap–
I mean red corrected uncertainty intervals. The trend is just the OLS trend.
The reason you are getting a different trends is I made a boo-boo related to GISTemp not being downloadable. The boo-boo is also related to programming inefficiency… that said, I think those graphs are NOAA/NCDC. See the update from yesterday afternoon.
nasaUrl <- "http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt"
download.file(url=nasaUrl,destfile=basename(nasaUrl),mode="wb")
works just fine for me on windows
mosher
download.file(url, destfile=file, mode="wb", method="curl")Worked. It didn’t work without the “curl”.
lucia (Comment #110562) your Marilyn Monroe trend with a 60 year peak/peak cycle is just what I got.
Doc–
That’s HaroldW’s graph. I just added the HTML so everyone can see the 30 year trends bugs wanted us to look at.
Lucia
file <- c(“GLB.Ts+dSST.txtâ€)
makes no sense.
c() is concatinate
so c("a","b","c")
gives you [1] "a" [2]"b" etc
you just want file <- 'GLB.Ts+dSST.txt'
or use the basename() function to strip off everything from the url
( works with directory structure as well ) to get the basename()
of a the path+file.name
you can also strip off the path, but that is another day, nother lesson.
RCurl as mentionedd above by mt is also an option. Try my fix. If that doesnt work, then I'll work up RCurl for you. It's far more powerful, but requires a bit more time to master.
mt might have his RCurl code handy, otherwise I'll have to dig into old code. The nice thing about RCurl is that you can suck down the directory structure and figure out where files are and get dates and sizes etc to that you can check for updates and check download size etc etc.
However, RCurl can be a PITA especially on MAC I have found. Its going to link to libcurl and so that needs to be up to date.
wget() has always been spotty for me on windows..so I stick with RCurl ( on curl obviously) for any fancy work with getting data.
Support from the maintainer ( I think its duncan?) is Ok, so questions should be posted on R list or maybe mail him directly
mosher–
I tried precisely your fix using what you wrote literally and it did not work. So, I looked up and tried a few things which did not work. This did work:
download.file(url, destfile=file, mode="wb", method="curl")Other things might work too.
I’m sure I did the concatenation of filenames by simply imitating someone else’s file. It may make no sense, but it seems to have nothing to do with the problem downloading from GISS. I tried “download.file(url, destfile=file, mode=”wb”, method=”curl”)” with or without the destfile=file, mode=”wb” and method=”curl. The first thing that worked is the bit I cut and pasted above.
HaroldW, I like the graph too but I’m wondering what you did to get the 30-year trend for 2012?
Carrick,
I think it is the trend for the 30 years leading up to the plotted year.
Steven #110573
I looked at RCurl, but they said there was no Windows version.
I couldn’t see that it offerred any advantage over just using curl on the command line. It’s much better documented. And it’s just a 350k executable.
To download a file only if there is a more recent version:
curl -z myfilename theirURL.
SteveF, that makes sense, especially if you look at the fit the 60-year + a nonlinear secular variation: Figure. (This was a fit I did to HadCRUT4 that we had discussed on a previous thread.)
Using the midpoint of the fitted intervals, I found the most recent maximum to be around 1990.
(But I’d suggest to HaroldW that modify his graph to show the mid-point!)
Carrick (Comment #110579)
Interesting graph and one I had not previously seen. How was the presentation smoothed?
LOL at justification for using shorter timeframes.
Can SteveF produce a 30 year trailing avg of the PDO or HaroldW produce a 25year OLS trend for the global warming rate? So we can see how the two graphs line up.
Please
(Before Lucia tells me to do it myself I’m going to claim technical incompetence)
nick i probably built it from sources. it just wraps curl.
the advantage is staying inside of R. that is to say there is no technical advantage. one can send email from within R. its an R nut thing. we like to do everything within R. principles nothing really technical
HR,
Here is the 30 year trailing average of the PDO index.
http://i53.tinypic.com/2hyexs7.png
Kenneth Fritsch, I just fit a fourth order polynomial + a sixty year cycle, then took the first derivative.
So fit:
$latex T(t) = T_0 + T_1 t + T_2 t^2 + T_3 t^3 + T_4 t^4 + a \cos \omega t + b \sin \omega t$
with
$latex \omega = 2\pi/ 60 \,\hbox{years}$,
then take the derivative…
I’m not claiming anything fancy here, but I did verify I wasn’t over-fitting the polynomial.
steven mosher:
It’d be interesting to compare download speeds…
Instead of trying to justify the various versions of station data, GISS would be more convincing if it simply proclaimed what it wanted each time series to be.
Carrick,
when lucia writes download.file() this is what happens.
basically download.file just wraps the console call (system) that Nick wants to make. When she doesnt specify a method and leaves it blank, it defaults to her internal method which is probably wget(). I’m guess her wget() is hosed on her machine, but her curl works great.
The way we look at source in R is just type in the function without the parenthesis.. and the source is dumped.
So its not like we rewrite curl or anything we just wrap it so you dont have to bounce in and out of environments.
> download.file
function (url, destfile, method, quiet = FALSE, mode = “w”, cacheOK = TRUE,
extra = getOption(“download.file.extra”))
{
destfile
method <- if (missing(method))
getOption("download.file.method", default = "auto")
else match.arg(method, c("auto", "internal", "wget", "curl",
"lynx"))
if (missing(mode) & length(grep("\\.(gz|bz2|xz|tgz|zip|rda|RData)$",
url)))
mode <- "wb"
if (method == "auto") {
if (capabilities("http/ftp"))
method <- "internal"
else if (length(grep("^file:", url))) {
method <- "internal"
url <- URLdecode(url)
}
else if (system("wget –help", invisible = TRUE) == 0L)
method <- "wget"
else if (system("curl –help", invisible = TRUE) == 0L)
method <- "curl"
else if (shell("lynx -help", invisible = TRUE) == 0L)
method <- "lynx"
else stop("no download method found")
}
if (method == "internal")
status <- .Internal(download(url, destfile, quiet, mode,
cacheOK))
else if (method == "wget") {
if (quiet)
extra <- c(extra, "–quiet")
if (!cacheOK)
extra <- c(extra, "–cache=off")
status <- system(paste("wget", paste(extra, collapse = " "),
shQuote(url), "-O", shQuote(path.expand(destfile))))
}
else if (method == "curl") {
if (quiet)
extra <- c(extra, "-s -S")
if (!cacheOK)
extra <- c(extra, "-H 'Pragma: no-cache'")
status <- system(paste("curl", paste(extra, collapse = " "),
shQuote(url), " -o", shQuote(path.expand(destfile))))
}
else if (method == "lynx")
status “, shQuote(path.expand(destfile))))
if (status > 0L)
warning(“download had nonzero exit status”)
invisible(status)
}
Nick.
I was mistaken I thought I built from sources or got it from omegahat, but looking back I got it from dr. Ripley.
http://en.wikipedia.org/wiki/Brian_D._Ripley
Since it depends on curl being installed they dont put it on CRAN
But Dr. R supports us by doing builds himself
http://cran.r-project.org/bin/windows/contrib/r-release/ReadMe
Packages
GLMMGibbs, RCurl, XMLRPC, survnnet, yags
do not build out of the box or other special in other circumstances.
Nevertheless these are available at
http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/2.15/
kindly provided by Professor Brian D. Ripley.
HaroldWHR-I told bugs to do it himself because I strongly suspect he was trying to “make a point by assigning work to others”. If he wants to make a point, it would be better to do so directly. Anyway, even he could presumably make a graph and fit a trend.
Carrick(#110575 & 110579) –
Apologies for being so late to respond, but when I try to post from work, it runs afoul of Lucia’s scripts…As SteveF suggested, the trends are plotted according to the end of the 30-year interval. So the points plotted at x=Dec2012 are the trends for the period Jan1983-Dec2012 inclusive.
.
I tried to indicate this by labeling the x-axis as “Year (end of interval)” but that’s in fairly small lettering and easily overlooked. While I like your suggestion of plotting the points against the midpoint of the interval (and often do just that), I was trying to avoid giving the impression of not being current, because the last point plotted would be at year=1998. If I update the chart, I’ll add a text box to the plot to make it more evident, whichever convention I use.
Lucia (#110592) –
That’s HR to whom you responded, not me.
@lucia (Comment #110564)
February 26th, 2013 at 9:38 am
Why?
It would be interesting to compare the ten year to 30 year trend.
Thanks, HaroldW, but what I don’t like about using the ending year is you can’t plot the 10-year curve (for example) on the same graph as the 30 year, plus IMO it’s misleading (it makes it look like you have results from a period where you don’t).
bugs–
You have the data for both now. Go ahead, compare them and reveal to the rest of us whatever you think is interesting about the comparison.
bugs: When something changes for a decade, a 30-year straight line wont reflect as much of a change as a 10-year straight line. I can predict that without comparing any data.
In fact, I’ve done simulations that showed me if you take temps from the 1950’s to date and add 100 years of unchanging temperatures (a little noise added to 2012’s temperatures, there’s still a statistically significant positive slope. Straight lines and time series work that way, and statistical significance is not practical significance.
Carrick (#110596) —
Agreed about the mis-registration which occurs if one does not plot a windowed average (or trend calculation) at its center. [I thought I could get away with it on that chart, as all curves have the same 15-year lag.] Your point is well taken, and here is the updated version.
.
Edit: Data still current only through Dec2012. Jan2013 values are not currently available for HadCRUT3&4.
Re: Carrick (Feb 26 20:09),
Why not fit the frequency too rather than specifying it? When I do that with the AMO index, I get a 67 year cycle.
DeWitt:
To be honest, I was just being lazy (as in, it didn’t matter for what i was doing)… but I have done both.
In any case, HadCRUT 4 is rather an odd duck (which you can see from HaroldW’s plot). I get 67 years for this series, but for NCDC & GISTEMP I get something around to 56 years. Ideally, you’d want to constraint the frequency of the oscillation using the fit to all three series using a robust estimator.
Doable, but more than this is worth with the limited data availability and the lack of a credible predictive. theory that can be tested with existing observations.
lucia (Comment #110597)
February 27th, 2013 at 6:40 am
I am not a statistician or mathematician.
You said I didn’t want to see something. What did you mean by that.
bugs,
I suspect you don’t want to see that we have a number of choices for the 30 year trend and that the 30 year trend ending now is lower than the 30 year trend ending several years back.
As for not being a statistician or mathematician: I should think if you were smart enough to graduate from the 6th grade, you could learn to use EXCEL and make a graph. If, after haunting climate blogs for years and claiming you would be interested in seeing graphs, you have not learned how to make such a graph, I suspect you are either:
1) Not smarter than a 5th grader (because otherwise you would teach yourself to make the graph)
2) or don’t really want to see the graph (because otherwise you would just make the graph you claim you want to see.) or
3) are really truly lazy and/or want to waste other people’s time asking them to make custom graphs you could perfectly well make yourself or
4) I am mystified.
In all cases, I don’t think it’s worth my time to make a graph for you. Because if (1) is true, you really aren’t smart enough to interpret graphs or know which graphs are worth making. If (2)-(3) I don’t see why I should waste my time making the graph. If (4)… well.. I’m just mystified, but I still don’t know why I should spend time making you the graph of your dreams. I suggest you make it yourself. Really!
download.file(url, destfile=file, mode=”wb”, method=”curl”)
This does not work for me. I get this error:
“Warning message:
In download.file(url = nasaUrl, destfile = basename(nasaUrl), mode = “wb”, :
download had nonzero exit status”
(I am using R 2.14 on Windows)
SRJ–
Interesting. Well… solved my problem… 🙂
I had an “issue” reading RSS this morning. I never had the problem before– the file seemed corrupted. It fixed itself. Sometimes I wonder about the web.
“It would be interesting to compare the ten year to 30 year trend.”
Actually, the climate plotter will do this. Here is the 30-year OLS trend for GISS over time, plotted (centred) in blue, and the ten year trend plotted in brown. It’s a kind of hack in the plotter – regress is offerred as a “Smooth”, which it isn’t, and changes the units. The axes are right in °C/yr, but the label is wrong.
Thanks Nick. I asked a climate scientist once about short term, (ten years or so), trends, and he said they were a waste of time. That graph seems to confirm that opinion.
I have no idea, though, how you did it with those parameters.
Also, the utility works pretty well on Chrome, but it didn’t run at all on IE. Is it supposed to work on IE?
Bugs,
Yes, it runs on IE but you may have to tell it to allow Javascript. To get the trend plots, first select the graph – eg GISS. Then go to the Smooth button, set the time (eg 30yr), and from the selection box, Regress. Then in the panel to the left, tick GISS, and click the Smooth button. You may have to adjust the scale of the plot using the vertical pink and blue bars. To make a second curve, select GISS again etc
Nick Stokes (#110635) interesting comparison. But if the data refer to the 30 or 10 years leading up to the plotted year,why are the values for recent years absent?
wilt–
People can plot based on ending year, beginning year or centered year. Centered year is most useful for comparing the 10 year and 30 year averages because the 10 year period is now in the “center” of the 30 year period. But when this is done, the most recent 10 year average is the one from Jan 2003-Dec 2012. It is shown as the centered average for the time midway between those to end points– so near 2007.5 which will be between june and july. The most recent 30 year average is further back in time.
Lucia, thanks for your clear explanation. I hadded missed the term “centered”.
SRJ do you have curl installed? As Steven Mosher explained R implements its curl routines via system calls.
lucia (Comment #110621)
February 27th, 2013 at 3:39 pm
You suspect wrong. The 30 year trend is a far better indicator of the warming signal than a ten year trend, as I have been told by an expert in the climate science area of knowledge. It looks like he was correct. Refer to the graph Nick was kind enough to plot for me. I will have to have a look at his utility more closely, it looks a lot more useful than Woodfortrees.
bugs–
Well, then you’ve gotten something you want to see. Good.
The purpose of the shorter trends is more to test projections. You can’t test projections based on data available prior to the projections being made and which contributed to the tools used to make the projections. So, for that, the 30 year trend is not particularly useful.
Lucia,
Sure you can! It seems to be done regularly by climate modelers… among many others. And damn if those ‘projections’ don’t usually end up looking pretty good.
.
Maybe your statement should be nailed to the door at GISS in Manhattan.
HaroldW
I hope you don’t mind but I took the liberty of (roughly) adding CO2 to your 30 year trend graph
https://lh3.googleusercontent.com/-XNCRVvsSGFU/UTIQqEiX80I/AAAAAAAABQo/fXVHHv9lFWU/s803/30yearTrends.jpg
Some questions which arise
1. Why did the temperature trend rise at all between 1910 & 1945?, during which time the CO2 trend rose from 0.2 ppm/yr to 0.25ppm/yr, with CO2 rising from 298ppm to 308ppm
2. Why did the 30 year temperature trend progressively fall to negative values between 1945 & 1965?, despite the CO2 trend accelerating to double the prior period
3. Why was the rate of rise of the 30 year trend in the period 1925-35 faster than at any equivalent decade in the modern warm period (1970-2000)?, despite the rate of rise in atmospheric CO2 in 1925-35 being trivial
4. Why has the 30 year trend post 2000 declined?, despite CO2 trend continuing to rise through the period
SteveF
I believe you may have jumped to a false conclusion, I agree that increasing atmospheric CO2 produces a positive forcing to global mean temperature. I also agree that global mean temperature has risen approximately 0.8degC since pre-industrial times.
I simply took issue with your statement that a positive CO2 forcing MUST result in a rise in GMT, as it is clear that CO2 forcing has been overwhelmed by other forcings/feedbacks in the recent past.
Gras Albert,
“I simply took issue with your statement that a positive CO2 forcing MUST result in a rise in GMT, as it is clear that CO2 forcing has been overwhelmed by other forcings/feedbacks in the recent past.”
.
So let’s see; are you suggesting that added forcing does or does not cause warming, all else being equal? If you are only suggesting that other things (like seasonal changes, or pseudo-cyclical patterns like ENSO) can have a greater influence short term, then sure. But that doesn’t really mean much. If you are suggesting that you are certain that rising GHG’s have not caused any of the observed warming, then you should be able to tell us mechanistically exactly how that could happen, considering that radiative physics indicates added GHG’s should cause some warming.
Lucia: A very late question, but you do a lot of work with trends. I thought that the whole ARMA(1,1) noise (or AR(1) + white noise) thing was to adjust confidence intervals and did not affect the actual slope at all. In fact, I’m pretty sure I’ve read the straight OLS yields exactly the same slope, it just underestimates the uncertainty.
As a first cut at ARMA(1,1) noise, I used R’s arima:
arima (x, c(1, 0, 1), xreg=time (x))
and I definitely get a different slope and intercept than using lm:
lm (y ~ x, data=data.frame (x=time (x), y=x))
Is arima not doing what I read it to do (with xreg)? Or am I mistaken about ARMA(1,1) noise not changing the slope? Or something else?
Wayne2
did not affect the actual slope at all
Depends on how you do it. If you fit ARIMA using R, the noise model affects the best estimate of the slope. If you widen error bars for ARMA11 “the Tamino way”, he uses the slope from OLS.
You’re mistaken about this. You probably read a description of the Quennoulle method or the Tamino ARMA method. Those just widen error bars.
We could argue about which method of estimating the trend is better (and to some extent ‘it depend’)– but you’ll notice the effect on the slope is small. In contrast, the effect on the size of the error bars is large. So.. sort of doesn’t matter.
Lucia, the question I have is does the ARMA(1,1) method change the expectation value of the slope?
Carrick–
If you use the ARIMA function to fit (y,x) instead of using Least Squares, the computed best fit trend is a bit different. I think the reason is that the ARIMA-fit function includes an attempt to find the initial condition for the noise based on the model. So… the trend changes a bit.
That ‘best fit’ value is supposed to be the best estimate of the expected value.
The true expected value is– of course– not really affected by any analytical process it simply exists. But or best estimate of the magnitude will change.
Does that answer?
Sort of …
I realize that the fit slope using ARIMA is going to be different for a given realization of the noise field than simple OLS.
My question is:
Is the mean value over a large number of realizations of the noise field (e.g. Monte Carlo) of the fit slope in these two cases the same?
If not, which is closer to the “true expected value”?
(You can test this question using a Monte Carlo approach of course.)
Carrick, lucia can correct me if I’m wrong, but my understanding is if you have the right ARIMA model, it will produce a more accurate estimate (on average) than a simple OLS regression. If you use a wrong model, it may be equally accurate or less accurate depending on how the model is wrong.
The reason is the ARIMA approach takes a step to reduce noise prior to fitting the linear regression. By reducing the noise, you get a more accurate fit if you have an appropriate model.
Of course, I’m self-taught on these matters so I could be way off!
Carrick–
Yes. If you generate data with a known trend synthetically and run a lot of test, both methods return the correct trend. So, both agree with each other and with “true”.
OLS will give you the smallest standard deviation about the “true” value you used to synthesize the data and is advantageous for that reason. I *think* the ARIMA fit gives you better estimates for the next few points– even though there is a larger standard deviation in the estimate of the trend. That’s why I said “it depends”
lucia, are you sure “both methods return the correct trend”? My understanding is autocorrelation can introduce bias into the OLS results. If so, the OLS regression will give less accurate results than an ARIMA fit with the right model.
I did a quick test of this in R, and it seems to confirm my understanding. I generated a sequence of autocorrelated noise, added it to a linear line then used both methods to generate a fit. I then repeated the process ten thousand times.
The ARIMA fit produced a more accurate result. The difference was minor as we’d hope with 10,000 series, but it was there.
Brandon
hmmm… I’m pretty sure I’ve gotten the opposite result. But… well.. maybe not. I’d have to rerun.
Lucia,
Based on your experience, are the OLS and ARIMA trends usually pretty close to each other, or can there be large discrepancies? I had always guessed that except for very short trends the primary influence of auto-correlation was an underestimate of the uncertainty range calculated using OLS, not a big error in the calculated slope from OLS.
SteveF–
They are usually pretty close. I could gin up some monte carlo and show cross-correlations– but I’d rather finish up other things. The big issue is the size of the uncertainty intervals.
lucia:
I can’t think of any reason an OLS regression should produce more accurate results. Having an appropriate noise model should not produce a less accurate result. I could perhaps see OLS being more precise (the results I have so far are inconclusive on that), but not more accurate. I might be missing something though. Here’s some janky code I tossed together to test this:
Let me know if I messed up something.
SteveF:
My experience is they’re pretty close to each other. In the brief testing I’ve done today, I haven’t seen a case where the difference was more than five percent. On average, the difference has been about half a percent. You might get large discrepancies with really short series or other unusual cases, but otherwise, they should be quite close.
Lucia: One more question, perhaps off-topic: have you tried working with per-month time series such as January SH GISS Temps 1880-2012, instead of monthly series? It seems to be much cleaner working with temps all from the same month (and hemisphere), and I did a quick calculation, but may have messed up my accumulations, with something like:
which (foreach (j=101:120, .combine="c") %do% sum (foreach (i=1:12) %do% 0 < abs (Reduce (`+`, sign (confint (lm (y ~ x, data=data.frame (y=giss.nh.df[j:133, i], x=time (giss.nh.df[j:133, i]))))[2,])))) == 0)Yeah, sometimes I end up with obfuscated code. But if I’m doing it correctly, it looks like you can go pretty far back in time before you get statistically-significant slopes for the trend of a particular month’s temperatures.
Wayne–
Terminolgy question: In words what is the difference between “per month” series vs. “monthly series.”
Please don’t substitute code for words. I have no idea what you intend to do above. If what you mean is find out whether the trend in all “january” temps is statistically significant…. why would you want to bother with this? What notion do you mean to communicate with the term “cleaner”? (Real questions btw.)
As far as I can tell, the main thing that happens when you use january only is you have fewer data points which reduces power. That’s never a good thing. But if you think you achieve something useful, say what it is. (And no… not by showing code. In words.)
Lucia: Sorry, not sure of the proper terminology. 12x the data (entire GISS series rather than just January’s) would allow 3.5x (sqrt (12)) the effect size — I think that’s the term — wouldn’t it? On the other hand, things like annual cycles are just distractions and looking at how January has trended over time, and how that compares to July over time, seems to be directly meaningful.
That’s what I mean by “cleaner”, looking at NH January’s trend and perhaps any features it shows and not dealing with yearly seasonality, etc.
I’m guessing that the reason I can go farther back in time and not find statistically significant slopes (for each month) is because of the reduced power from having 1/12 the points, and not some real effect. Oh well.
Lucia: Another thought about “cleaner”. When I look at a plot of the partial autocorrelation for GISS NH temps 1880-2012, there’s significant correlation for the first 12 months. (I say “significant” meaning that the bars exceed the blue-dashed white-noise-based confidence interval as plotted by R’s ‘pacf’ function.) Would that indicate that we’re getting a lot less power out of the entire dataset than having 12x the amount of data (i.e. all months) would normally give?
Wayne–
I don’t want another thought about “cleaner”. I want to know what “cleaner” is even supposed to mean in the first place.
No. You can’t check for power by data mining for significance and finding the data set that triggers a “significant” over the minimum amount of time.
You check for power by running monte carlo simulations of the method. There may be advantages to throwing away data willy-nilly. But you are never going to achieve greater power by finding excuses throwing away good data. It can be ok to throw away bad, inaccurate, data. For example, if you did an experiment and noticed a thermocouple broke, that’s bad data and should be tossed. But “temperature from Feb-Dec” is not the criteria for “bad” data.)
There might be better ways to handle Jan-Dec data than deseasonalizing the anomalies and fitting the trend. But the better way is not to throw away Feb-Dec and look at trends for January only.
Opps… sorry. I missed your definition of cleaner!
Anomalies are deseasonalized to eliminate (or at least reduce) the annual cycle. The only real issue is whether it might be warming more in one month than another. If you want to figure that out, it can be worth looking at each month individually. But if you are trying to see overall warming, I don’t see any obvious advantage to looking at each month individually. Someone might come up with a method. If so, it would be interesting to see. But, since I don’t anticipate any big advantage to such a method, I’m not planning to try to develop such a method nor to hunt it down.
Yep. That’s why I don’t see this method as very promising. You might be able to gin something up based on using all 12 fits… but I’ll leave that work to someone who thinks it’s promising. 🙂
Lucia: Thanks so much for your responses. I’ve been curious about these things and don’t put much stock in whoever might answer a question in an open thread on WUWT, etc. But I do trust your insights.
I’m not saying that January NH GISS alone is better than all of NH GISS. Perhaps my idea is equivalent to throwing data out, but I view it more as treating GISS NH as panel data rather than as a single, univariate series.
I’ll look at panel-data-style analysis techniques to see if anything seems applicable. (Methods that come to mind, but may be useless include things like Hierarchical linear models, VAR, and GEE.) And also at running MC-style simulations to try to illuminate if working with months makes any sense.
Wayne2–
Hate to ask this.. but ‘applicable to answering what question‘? Or maybe I should say “to testing what hypothesis”?
Whatever it is you are trying to learn, setting up alternate methods of testing and running MC is a good idea. You can learn the power that way– if you test properly.