September GISS Temp up by 0.13°C

GISS LOTI went from 0.61°C in August to 0.74°C in September. because of the US shutdown, the report was exceptionally late, and in fact the surface indices were reported in the reverse of the usual order, with HADCRUt first. GISS showed the highest rise; NOAA 0.03°C, HADCRUT 4 and TempLS less than 0.01°C. But the satellite indices showed rises comparable to GISS.


Here is the GISS map for September 2013:



And here, with the same scale and color scheme, is the earlier TempLS map for September:

And here is the WebGL map of actual station temps during the month.

Previous Months

August
July
June
May
April
March
February
January
December 2012
November
October
September
August
July
June
May
April
March
February
January
December 2011
November
October
September
August 2011

More data and plots

New page - climate blog index

I described here an experiment I tried in July (when Google Reader faded). I selected about a dozen popular climate blogs, for which I regularly looked at the RSS files (every two hours). They aren't all blogs that I agree with, but they are ones at which I sometimes comment.

A database of posts and comments has accumulated, and I've done some organising to make that accessible without a huge initial download. So I've scrubbed it up and promoted it to a page, visible upper right.

It's just a table of blog names, thread titles, dates and authors (comments and posts). The dates are linked to the source. You can order by each column, or you can subset by each of those properties (in combination, if you want). You can, for example, put up a linked list of all your own comments on those sites over a time period.

It starts with just the current and previous month. But you can vary the time range (more data =slower). More details are at the page.

TempLS global temp up 0.01°C in September

GHCN now seems to be OK after the shutdown. There was a tiny increase from August (0.49°C) to September (0.504°C). UAH seems to be the only other index reporting to date - it showed a much larger increase.

Here is the spherical harmonics plot of the temperature distribution:



Warm spots in Canada and Australia; cold in Siberia and E Europe..


And here is the map of stations reporting:

Shutdown - TempLS delayed

GHCN has been updating, as has ERSST. But when I look into GHCN, oddly, there are only US data, no ROW. So no TempLS post yet.

FWIW, the data that is there suggests a slight fall in Sept.

The Shutdown and updating

Just a note on how the US government shutdown affects updated data here. Hopefully all will be fixed when it ends.

• NOAA SST data isn't getting through, so the SST movies and WebGL trackball daily plots have no October data.
  
• The NASA and NOAA websites are down, so there probably won't be a monthly index update for the duration. UAH is out and RSS should be OK.
  
• GHCN data is getting through, so TempLS should be OK - may be the only one until HADCRUT comes out. The monthly updates of surface temp should be OK.
  
• JAXA is unaffected.
  

Active ocean acidification calculator.

I wrote a post a few weeks ago on ocean acidification. The message was that there is more to acidification than pH, and that it's basically an issue of equilibrium between carbonate species. Difficulties of measuring pH can be avoided.

In the course of discussion, I said I would try to develop an active calculator to display some of the principles.

It calculates equilibrium adjustments, and I expect that is a rather specialised interest. I'm putting it forward to demonstrate

• that any two of the species are sufficient to determine the others - pH is not essential and can be treated as just a following variable
• CO2 when dissolved is almost completely reacted (at usual pH)

The main process of interest is the addition of CO2, which causes carbonate to go down. This affects the solubility of CaCO3. I haven't shown that equilibrium, because it only makes sense where solid CaCO3 is present. The calculator has two modes. You can show the effect of adding a fixed amount of any reagent, or you can specify any two. Details and gadget below the jump.

The gadget

There are four species (H+,CO2,HCO3- and CO3--). For details of the chemistry, see my previous post. I also have two combinations, dissolved inorganic carbon (DIC) and total alkalinity (TA). These are conserved during the reaction as total C and nett charge. There are six bars, four in gold for the species, and two for the combinations. They are scaled on a negative log scale (like pH). You can click on each bar to vary each constituent. But first you need to set a constraint using the radio buttons at the top. The default setting is marked "Add". Here you can only vary the species (gold). Your click will indicate what the concentration would have been if the change (+ or -) happened without reaction. But in fact your addition is added appropriately to DIC and TA, and then the equilibrium adjusts with these held constant. Sometimes almost all what you have added reacts, and you see only a small change in that variable. If you click another radio button, you move to "set" mode. The column you have marked is held fixed, and the column you click (a different one) will move to the value you have set (if possible - see below). "Set" is a slight misnomer - it's better seen as the values the other components are known to have when the "set" values have been measured. The blue bars show the current state, with log values beside. The absolute values show in the table on the right, with the pK values. The intent, not yet implemented, is that you can vary these.

Update: You can now vary pK1 and pK2 (but see Notes below). Modify the text box and click the new "Accept Data". You can use this to look at totally different equilibria if you want. To make it monoprotic, just set pK2 to, say, 15.

Numerics

I'm using a Newton-Raphson process to solve the equations. This is fast and accurate, but may not necessarily converge. On the right there is a purple box which tells you whether your request has converged or failed. If failed, for whatever reason, no change has been made and you can try again. Some requests are impossible. For example, no carbon species can exceed DIC. But sometimes the Newton process just had trouble with the size of change you asked for. Try getting there in two steps.

Graph

A Bjerrum plot is shown (from the previous post). It shows equilibrium concentrations normalised against DIC vs pH. A thin gold bar moves to show where you currently are on the plot.

Here it is:

Notes

The constants pK1 and pK2 come from Zeebe's review article referred to in my previous post. I do not distinguish between CO2 and H2CO3, and the constants reflect that. [CO2] means the combined value. A pK1 based on carbonic acid alone would be lower, but you would then need the equilibrium between CO2 and H2CO3.

Adjusting temperature series stats for autocorrelation.

I've been writing a lot lately on autocorrelation for temperature time series. This is basically a methods post.

I've often referred to a post by Hu McCulloch, in which he reproves a scientist for quoting just a OLS regression for a monthly series with OLS CI's, and shows how the CI's are much larger if the residuals are assumed to be not white noise but AR(1) correlated. He gives a useful set of formulae, and describes in particular the Quenouille AR(1) correction to the trend uncertainty. This seems to be the current standard, and it is what I use in my trend viewer.

Stats program R has an arima(xreg=) function which will give trends with error estimate for general ARIMA models. It solves a slightly different optimisation for regression, with a different trend for each ARIMA type. But for AR(1), the results are similar to the Quenouille approach.

So why bother with approximations? I'm interested for two reasons, in connection with the trend viewer. In my weekly updates, I do about 5 million regressions. That takes a bit over an hour, but only because I avoid the individual regression summations, instead using differences of cumulative sums calculated for each dataset. You can't to that in a package. arima() would take many hours.

The other reason is that to provide interactive output, I have to be able to repeat the calculation in Javascript. Time isn't an issue there, but programming simplicity is.

In this post, I'll describe some extensions which allow Quenouille style corrections to more complex ARMA() models. I also give corrections to the OLS trend. The trend corrections are O(1/N) (N=#data) and I believe the error is then O(1/N^2), and for the CI's, O(1/N). Those errors are for the approximation to the stochastic model - whether it is the right model is a greater uncertainty.

With these formulae the complete procedure is

1. Do a normal OLS regression.
   
2. Calculate the residuals and estimate the first few autocorrelations (one lag for A(1)m 2 for AR(2) and ARMA(1,1) etc)
   
3. Use the Yule-Walker equations to estimate the ARMA parameters. That just solves a few linear equations in the AR() parameters.
   
4. Correct the trend and uncertainty using thge formulae given here.
   

Fitting stochastic models

I'm talking about fitting to regression residuals, so it is simpler to have a regression already done, although I described a method for simultaneously fitting regression and residual model. An ARMA model for regression residuals can be written: A(L)d=M(L)ε Here A is a polynomial in the lag operator L, A=1-A1*L-A2*L^2... applied to the time series d (so L d_i = d_{i-1}) M(L) is a similar polynomial in L: M(L)=1+M1*L+... ε is a vector of iid random variables The series d is taken to be the residuals of a regression on data y: d=y-a-b*x. The notation ARMA(j,k) means that A is of order j and M of order k. The effect of applying the polynomial A to d can also be represented by the lower triangular matrix (AR(2) example) A=

...0	-A2	-A1	1	0	0...
...0	0	-A2	-A1	1	0...
...0	0	0	-A2	-A1	1...

This is of course a fragment - if d has length n then A is (n-j) x n and is Toeplitz. M also corresponds to such a matrix. The regression equation can be written M^-1L d = ε Comprehensive programs like the R arima() will minimise d^* K^* K d where K = M^-1L over coefficients a and b, and return the appropriate uncertainties. The alternative described above is to minimise d^* d, and then adjust the coefficients and their uncertainty.

Finding and adjusting trend

The regression of the ARMA model gives trend b = y^* KK x/ x^* KK x, (KK = K^* K) assuming K x has zero mean. But KK is symmetric, and also near Toeplitz, deviating near the top left and bottom right corners. The i'th row of KK (far from corners) is symmetric about the i'th value (with unsymmetric truncation). But the arithmetic progression x can be split into a constant x_i and a part x-x_i that is antisymmetric about i. The latter yields a zero sum with that row of KK, so the result is R*x_i, where R is the (constant near centre) rowsum of KK. So KK * x ≈ R * x This argument gives b = y^* KK x/ x^* KK x ≈ y^* x/ x^* x which is the OLS value. That's why only a small adjustment is needed - it is due to end effects. In the case of AR(n) models, the deviation of KK from Toeplitz only affects the top n and bottom n rows, so it is a fairly simple matter to add up the missing terms. The corrections to the OLS trend for the various models are, for time (=x) series y of length N:

OLS or AR(0)	b
AR(1)	b+(3*b+6*(y_N*A1+y_1)/(1-A1)/(x_N-x_1))/N
AR(2)	b+6*(b+(y_1*(1-A1)+y_2+A2*y_(N-1)+y_N*(A1+A2))/(1-A1-A2)/(x_N-x_1))/N

And here are the results comparing these approximations with values computed with the R arima() program, setting xreg and using method="CSS".The corrected values are much closer than the OLS trend. Trends are from 1980 to July 2013.

Data	Hadcrut4	GISSLO	NOAA	UAH	MSU.RSS
OLS	1.5643	1.5626	1.4878	1.4101	1.2979
AR(1)	1.5743	1.5661	1.4936	1.4257	1.307
AR(1) approx	1.5742	1.5662	1.4936	1.4257	1.3073
AR(2)	1.5946	1.5909	1.5072	1.4705	1.3392
AR(2) approx	1.5959	1.5907	1.5079	1.4696	1.3389

Trend uncertainties.

The original adjustment for AR(1) autocorrelated regression residuals was given by Quenouille (1952), and has been widely used. In an earlier post, I gave the corresponding approximations for ARMA(1,1) and AR(2). I'll exrend to AR(3) and give the derivations here. Note that whereas for the trend we were looking for an O(1/N) adjustment, here it is O(1) - ie much bigger.

So it is sufficient to find an uncertainty adjustment for either the OLS trend or the true AR() trend. OLS is simpler.

Our ARMA model is
A(L)d=M(L)ε ε is iid
or, reverting to matrices
d=A^-1M ε
The variance of the OLS trend b is then estinated by
var(b)= x^*C x/x^* x
where x is time (mean removed) and C is the estimated covariance of d. In the model, C=M^*A^*-1A^-1M.
Now C is symmetric, and apart from end effects, Toeplitz, reflecting stationarity. That means that each row is symmetric about the diagonal. So when the i'th row of C multiplies x, the effect of the part that is odd about x_i makes no contribution. But the even part is just the constant . So the effect of applying C to x is just to multiply by the constant rowsum of C (apart from O(1/N) end effects). The ignoring of end effects here requires that there is only moderate autocorrelation, so C is diagonally dominant.

Each row of C is just the diagonal value (zero-lag variance of d) multiplied by the acf. That diagonal value divided by x^*x. So the correction required to the variance of b is just (with our 1/N caveats) the sum of the autocovariance function.

Now it is good to go back to polynomials. The (Toeplitz) matrix product M^*M is equivalent to the polynomial M(L^-1)*M(L). And the acf is almost just the set of coefficients of the Laurent series got by expanding M(L^-1)*M(L)/(A(L^-1)*A(L)) in powers of L. A Laurent series is just a power series with both positive and negative powers.

I say almost, because the zero'th order coefficient of that series is not 1. So we have to find it and divide by it. This could be done with messy series expansion calcs, but it is neater with complex analysis, replacing L by z.

The Laurent series may converge in an annulus in the complex plane, and here it converges on the unit circle. To get the central coeficient, divide by z and integrate on that contour. It's a rational function, so that just means summing the residues inside the unit circle.

How to do that? For the previous post I did the algebra by hand, but for AR(3) I resorted to Mathematica, and I'll describe that process. I write A(z)=(1-az)(1-bz)(1-cz), so
1/A(1/z)/A(z)/z=z²/((z-a)(z-b)(z-c)(1-az)(1-bz)(1-cz))
I don't need to ever find the roots a,b,c, and they don't have to be real. To get the residue at a, I define f(a,b,c)=a²/((a-b)(a-c)(1-a^2)(1-ba)(1-ca))
Then I ask Mathematica to Simpliify[f[a,b,c]+f[b,c,a]+f[c,a,b]]
It comes back with a polynomial in the numerator and factors in the denominator, symmetric of course in a,b,c. The factors in the denominator I can reassemble as values of A(z) at specific points. The numerator I can recognise symmetrics like ab+bc+ca, and re-express these as coefficients of A.
Let A0 be the result - the zero'th order coefficient of 1/A(1/z)/A(z). Then the numbers in the acf are just the Laurent coefficients of 1/A(1/z)/A(z)/A0, and the sum - the Quenouille correction - is just 1/A(1)/A(1)/A0.

So here are some results. First the formulae for the Quenouille factors for adjustment of trend uncertainty:

AR(1)	(1+A1)/(1-A1)=A(-1)/A(1)
ARMA(1,1)	A(-1)/A(1)*(1+M1)2/(1+M12+2*M1*A1))
AR(2)	A(-1)/A(1)*(1+A2)/(1-A2)))
AR(3)	A(-1)/A(1)*(1+A2+(A1-A3)*A32)/((1-A2-(A1+A3)*A32))

And here are the numerical results for the monthly data, period 1980-July 2013. First the OLS trend and standard error; then the various models, with in each case the R arima() trend and se, and then the se as calculated by the above Quenouille-style approx.

Result	HADCRUT4	GISS	NOAA	UAH	MSU-RSS
OLS trend β	1.7020	1.7092	1.7500	1.4086	1.5676
OLS s.e. β	0.0463	0.0651	0.0514	0.0665	0.0641
AR(1) trend β	1.6984	1.7084	1.7449	1.4161	1.5678
AR(1) s.e. β	0.0785	0.0936	0.0765	0.1101	0.0978
AR(1) Quen s.e.	0.0781	0.0933	0.0762	0.1096	0.0974
ARMA(1,1) trend β	1.6953	1.7025	1.7423	1.4188	1.5687
ARMA(1,1) s.e. β	0.0946	0.1132	0.0936	0.1189	0.1027
ARMA(1,1) Quen s.e.	0.0940	0.1123	0.0930	0.1182	0.1022
AR(2) trend β	1.6854	1.6831	1.7363	1.4221	1.5756
AR(2) s.e. β	0.0929	0.1103	0.0899	0.1228	0.1045
AR(2) Quen s.e.	0.0925	0.1106	0.0894	0.1217	0.1037
AR(3) trend β	1.6893	1.6853	1.7395	1.4098	1.5673
AR(3) s.e. β	0.0908	0.1060	0.0901	0.1073	0.0971
AR(3) Quen s.e.	0.0900	0.1059	0.0893	0.1060	0.0961

Conclusion

For autocorrelated data, the ordinary OLS trend analysis will get the trend near right with a O(1/N) error, which can be calculated. But the standard error will be substantially underestimated. There are at least three options

• Use a package like arima() in R, which probably uses a general optimisation algorithm
  
• Use the Newton-Raphson-style algorithm that I described here (I'm planning an update with more systematic code). This is fast and does everything at once. It is easy to code in R, but not so much in, say, Javascript.
• Calculate and correct the OLS trends, using Yule-Walker equations to estimate model coefficients.

The approximations will deteriorate with vary high autocorrelation. Then the best remedy is to instead use, say, annual averages. Little information will be lost, and the autocorrelation will reduce.