Climate Scientists and the Riemann Hypothesis

Jul 7, 2008 – 6:43 AM

During the last few days, there has been a flurry of activity following an announcement by Xian-Jin Li of a supposed proof of the Riemann Hypothesis, perhaps the most famous unsolved mathematics problem (now that Fermat’s Last Theorem has been solved.) Atle Selberg, an eminent mathematician, recently said:

If anything at all in our universe is correct, it has to be the Riemann Hypothesis, if for no other reasons, so for purely esthetical reasons.

On July 1, 2008, Li posted his proposed proof on arxiv.org here, stating in a short abstract:

By using Fourier analysis on number fields, we prove in this paper E. Bombieri’s refinement of A. Weil’s positivity condition, which implies the Riemann hypothesis for the Riemann zeta function in the spirit of A. Connes’ approach to the Riemann hypothesis.

I recently noticed a quip that there were two sorts of mathematics papers – those where you can’t understand anything after the first page and those where you can’t understand anything after the first sentence. I’m glad that I’m not alone.

Subsequent events showed an interesting role for blogs even in a mathematical field where there are only a few people in the world who can read the paper.

Terence Tao observed the next day (July 2) on his blog (as a comment not even a thread):

It unfortunately seems that the decomposition claimed in equation (6.9) on page 20 of that paper is, in fact, impossible; it would endow the function h (which is holding the arithmetical information about the primes) with an extremely strong dilation symmetry which it does not actually obey. It seems that the author was relying on this symmetry to make the adelic Fourier transform far more powerful than it really ought to be for this problem.

During the next couple of days, Li put two revised versions of his preprint onto archiv.org and seemed to feel that he could cope with this problem.

On July 3 at another blog, Fields Medalist Alain Connes stated:

I dont like to be too negative in my comments. Li’s paper is an attempt to prove a variant of the global trace formula of my paper in Selecta. The “proof” is that of Theorem 7.3 page 29 in Li’s paper, but I stopped reading it when I saw that he is extending the test function h from ideles to adeles by 0 outside ideles and then using Fourier transform (see page 31). This cannot work and ideles form a set of measure 0 inside adeles (unlike what happens when one only deals with finitely many places).

On July 6, it was announced on arxiv.org that the paper had been withdrawn.

This paper has been withdrawn by the author, due to a mistake on page 29.

Li’s failed proof was 40 pages of dense mathematics, with one flaw on page 29. But the flaw was enough to cause the paper to fail.

Just think how much easier it would have been for Li had he published in Nature. His article would have been 3 pages long – a little section on the colorful history of the Riemannn Hypothesis, moving quickly to the “results” and implications. Maybe Nature editors would suggest that he mention the number of zeros was unprecedented. Due to “space limitations”, there would obviously be only a few sentences in the running text on how his proof actually worked, but the words “rigorous” and/or “conservative” would almost certainly have been used.

Their running text in Nature would perhaps said:

“We used the E. Bombieri’s highly conservative refinement of A. Weil’s rigorous positivity condition, which implies the Riemann hypothesis for the Riemann zeta function.”

Isn’t that much more convincing than Li’s original? Now you can see the realdefect in Li’s proofs: the failure to use either of the terms of “rigorous” or “conservative”. Instant climate science rejection. Had he correctly used the terms “rigorous” and “highly conservative”, nothing else would need to have been said and the paper would, of course, become an instant classic in climate science. Nobody would actually read it past the abstract, but it would be highly cited.

Had he published in Nature, if asked to provide details on his proof, Li could have refused and been backed up by Nature’s editors, who would have said that their policy does not require authors to show minute details of their methodology. If anyone doubted Li’s proof, if they were so smart, why didn’t they prove the Riemann Hypothesis themselves?

Perhaps 6 years later, someone would decode entrails of the Li “proof” from snippets, whereupon Li would haughtily announce that he had “moved on” to a new proof using an even more rigorous, conservative (and opaque) methodology and so he had been right all along.

By Stephen McIntyre | Posted in General | Comments (201)

Sea Ice Utilities

Jul 6, 2008 – 5:14 PM

I’ll post up some utilities for reading sea ice here (soon). For now, I want to document some relevant aspects of the data sets.

NSIDC Daily Data Sets
Arctic sea ice information from NSICD is available in two directories:
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north – has daily 2008 data;
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/final-gsfc/north/daily – has daily data for years prior to 2008 with each year having its own subdirectory

Satellites
In order to download a daily record, you need to be able to specify the satellite id, as these are incorporated in the directory names. Since 1979, there have been 5 satellites: n07, f08, f11, f13 and f15. During most periods, only one satellite is operating, but there have been overlaps of at least a couple of months. In 2008, there are two satellites: f13 and f15, with f15 coming into service only this year. Here’s a summary of relevant information (jstart, jend being julian days since 1970-01-01.)

sat start end jstart jend
n07 1979-01-01 1987-08-20 3287 6440
f08 1987-07-09 1991-12-31 6398 8034
f11 1991-12-03 1995-09-30 8006 9403
f13 1995-05-03 2008-12-31 9253 14244
f15 2008-01-01 2008-12-31 13879 14244

I’ve made a function read.arctic (day,month,year,sat=”override”) to read the daily binary data. The override satellite specificaiton is the most recent satellite if there’s a choice, with the option of specifying the older satellite. This function returns a 304×448 array of sea ice concentrations (with a few special codes).

Cell Areas
I’ve been unable to locate a file showing the areas of the gridcells (which are said to vary) and with the variation affecting the calculation of sea ice area. Daily sea ice area appears to be calculated as the sum of the area of gridcells with sea ice concentrations greater than 15%.

Information on the gridcell structure is provided here . On the left is a diagram from that reference; on the right is a plot of downloaded sea ice information ( downloaded in a 304×448 matrix, with the 448 column then reversed so that the plot function matches.) The coordinates of the corners are thus available from this diagram. Note that the rectangle is slightly asymmetric and that the north pole is not at the center of the rectangle (though it’s close).

 polar2.gif  polar1.gif

The form of the plot on the left appears to be a “polar stereograph”. Whether I’ve named it correctly or not, each 10-degree circles seems to be uniformly spaced, so the projection does not preserve area. The url mentioned above says:

The polar stereographic projection specifies a projection plane or grid tangent to the Earth at 70 degrees. The planar grid is designed so that the grid cells at 70 degrees latitude are 25 km by 25 km. For more information on this topic please refer to Pearson (1990) and Snyder (1987). This projection often assumes that the plane is tangent to the Earth at the pole. Thus, since there is a one-to-one mapping between the Earth’s surface and the grid at the pole, there is no distortion at the pole. Distortion in the grid increases as the latitude decreases, because more of the Earth’s surface falls into any given grid cell, which can be quite significant at the edge of the northern SSM/I grid where distortion reaches 31 percent. …. To minimize the distortion, the projection is true at 70 degrees rather than at the poles. This increases the distortion at the poles by three percent and decreases the distortion at the grid boundaries by the same amount. The latitude of 70 degrees was selected so that little or no distortion would occur in the marginal ice zone.

It doesn’t seem like it would be that hard to construct an equiangular projection. According to my interpretation of this, the relative distortion from the projection used here is given by:

K= (90-70)*pi/180 / (sin ((90-70)*pi/180)) ;K
f=function(lat) K* sin ((90-lat)*pi/180) / ( (90-lat)*pi/180 )
#f(c(85,70,60,45)) #] 1.0193054 1.0000000 0.9746015 0.9188631

According to my calculations, the distortion at the North Pole is 2.06% and at the most remote of the grid is about 15%. I can’t replicate their distortion numbers, but the calculation is only high school geometry and my calculation looks sound to me.

I’ve made a function to calculate sea ice area incorporating the above information. So, here’s a pretty way in which you too can calculate your own daily sea ice area for a given day.

Calculating Daily Sea Ice
Here’s a nice simple way that you can calculate daily sea ice area for yourself:

source(“http://data.climateaudit.org/scripts/seaice/nsidc.functions.txt”)
day=5;month=7;year=2008
test=read.arctic(day,month,year);dim(test)
area_seaice(test) # 9.57941

Here’s a way to plot this (requires the package fields, available from http://www.r-project.org):
library(fields)

plot.arctic(test)
title(paste(“Sea Ice: 2008″,month,day,sep=”-“))

The satellite may affect measurements – a factor that I haven’t seen discussed at length, though obviously it’s the sort of thing that the specialists keep track of and do the best they can. In 2008, both the f13 and f15 satellites are up. Here’s the corresponding f13 area:

test.f13=read.arctic(day=5,month=7,year=2008,sat=”f13″);
area_seaice(test.f13) # 9.41061

Since 2007 measurements were done with f13, which measured about 160,000 sq km less than f15, a consistent f13 comparison would close the gap commensurately. These numbers are close to, but don’t reconcile precisely to published daily numbers. I haven’t located exactly how the daily numbers are calculated. In my “Race Reports”, I’ve been using daily totals from JAXA, which don’t seem to tie in exactly to NSIDC.

By Stephen McIntyre | Posted in General | Comments (22)

New Ice and Old Ice

Jul 4, 2008 – 4:24 PM

I’ve got a question for readers: can someone explain to me how the NSIDC sea ice projections were actually made? I can’t get from A to B. Relevant information is in their monthly reports for April 2008 and May 2008.

Remember how much we’ve heard about new ice and why this was expected to lead to 2008 ice meltback surpassing the 2007 record. Let’s try to follow their projection method step by step,

First, we start with March ice, which is the yearly maximum. They say:

Average sea ice extent for March 2008 was 15.2 million square kilometers.

Next we learn:

Figures 4 and 5 indicate that relatively thin, first-year ice now covers 72% of the Arctic Basin, including the region around the North Pole; in 2007, that number was 59%.

By simple multiplication, the area of new ice in March was 10.94 million km^2 (old ice: 4.26 million km^2).

Next we learn:

about 30% of first-year ice typically survives the summer melt season, while 75% of the older ice survives

Their projection using average survival rates is reported as follows:

If we apply the survival rates averaged over all years to current conditions, the end-of-summer extent would be 3.59 million square kilometers.

Can anyone explain this to me? If I multiply opening new ice (10.94 MM sq km) by 30%, I get 3.28 MM sq km; and multiplying opening old ice (4.26 MM sq km) by 75%, I get another 3.19 MM sq km, for a total of 6.48 million sq km – a figure which, by the way, looks like a pretty good estimate to me of where we might end up this year. If you inadvertently divide the surviving old ice by 10, you get something pretty close to 3.59 million sq km ( 3.28 + 3.19/10).

I don’t pretend to be familiar with these calculations and I could easily have been wrongfooted somewhere along the way. So if anyone can figure the calculation out, I’d appreciate it. BTW wouldn’t it be nice to be able to simply look at an archived script and be able to see for yourself what they did?

By Stephen McIntyre | Posted in General | Comments (53)

Climate Audit and NOAA FOI Policy

Jul 3, 2008 – 9:38 AM

For some reason, Michael Tobis seems to think that the Team is busy cranking out responses to an endless stream of my data requests. Nothing could be further from the truth. For me, getting data from the Team is no easier than posting at realclimate. At this point, it seems to be Team policy to simply not provide any data to me. So it doesn’t take them any time at all.

I racked my brains trying to think of whether I got any data from the Team last year (and by Team here, I mean the core Hockey Team, not the IPCC). The only thing that I could think of was the list of stations used in Jones et al 1990, which was eventually provided under a UK FOI request. I encountered an interesting little back story about this request today. Continue reading →

By Stephen McIntyre | Posted in Archiving, FOIA | Tagged , | Comments (39)

E-Mail, “Personal” Records and Privacy

Jul 2, 2008 – 5:32 PM

Recently, David Holland reported the unedifying spectacle of John Mitchell, Chief Scientist of the Hadley Center, an institution which proclaimed itself to be the “most significant” contributor to WG1, attempting to circumvent mandated requirements that IPCC be “open and transparent” and that all “written comments” be archived, by claiming that his email correspondence with IPCC Authors using his Hadley Center email account were somehow his “personal” property and therefore not subject to UK Freedom of Information legislation.

Elsewhere, we’ve seen the equally unedifying spectacle of Caspar Ammann, a U.S. federal employee, an employee of a federally funded research and development center (NCAR), seeking to circumvent the IPCC requirement that comments be “open and transparent” by sending comments to IPCC author Briffa at his CRU email, rather than through the IPCC registry (as even Susan Solomon did). Even though IPCC is required to archive all written comments, Ammann and Briffa regard themselves as above IPCC policy and entitled to exchange confidential comments, sort of like little school children whispering at the back of the classroom.

While the spectacle is obviously unedifying in itself, I thought that it would be interesting to examine government policies on e-mail correspondence to see whether there is any policy basis for either theory. Continue reading →

By Stephen McIntyre | Posted in climategate, FOIA | Tagged , , | Comments (29)

June 2008 Satellite Results

Jul 2, 2008 – 2:07 PM

You have to get up pretty early to be first out of the blocks on monthly temperatures. This month, Climate Audit is first out of the blocks with June 2008 monthly temperatures.

June 2008 MSU results (anomaly deg C), coming soon at http://vortex.nsstc.uah.edu/data/msu/t2lt/tltglhmam_5.2, are GLB: -0.114 (June 1988: 0.100), NH: 0.004 (June 1988: 0.140), SH: -0.232 (June 1988: 0.060), Tropics -0.334 (June 1988: -0.090). Here is a plot with the current month and Hansen’s June 1988 testimony month shown as red dots.

june2037.gif

Update (evening July 2): RSS June 2008
June 2008 from RSS shows: 70S-80N 0.035 (June 1988 0.091); SH -0.089 (June 1988 0.012); NH 0.154 (June 1988 0.166); Tropics: -0.116 (June 1988 -0.141). So UAH shows the tropics as going down by 0.24 deg C (RSS by 0.02) and similar large discrepancies elsewhere.

june201.gif

By Stephen McIntyre | Posted in General | Comments (120)

Sea Ice – the Stretch Run

Jul 2, 2008 – 8:03 AM

For anyone who’s betting that 2008 meltback will exceed 2007 meltback, I think that you’ll be able to pretty much know where you stand by the end of this week and your chances are not looking good right now based on this week’s exit polls. Another Climate Audit first. Continue reading →

By Stephen McIntyre | Posted in General | Comments (427)

Conrad Black and Mammoth Lakes CA

Jul 1, 2008 – 6:39 PM

The fraud trial of Conrad Black has been a very large story in Toronto. Black, a business tycoon, was recently sentenced to 78 months in Chicago; last week his appeal was rejected – news story here; blog account here; judgement here. Interestingly, the judge has his own blog.

Black and his co-defendants were executives of Hollinger Inc., a newspaper conglomerate (that, as a company has a long and interesting history, originating with the Hollinger gold mine in Timmins, historically the largest gold producer in North America.) As executives of Hollinger, they caused a subsidiary (APC) to pay non-compete fees of $5.5 million to themselves personally. APC had disposed of various newspapers leaving itself with cash and only one newspaper – in Mammoth Lakes, California, a town of 7903 people in the Sierra Nevadas, a town possibly familiar to Anthony Watts and some other CA readers in northern California, who will doubtless be surprised to learn of its prominence in such an important fraud trial.

From the judgment:

When it had only one left—a weekly community newspaper in Mammoth Lake, California (population 7,093 in 2000, the year before the fraud)—defendant Kipnis, Hollinger’s general counsel, prepared and signed on behalf of APC an agreement to pay the other defendants, plus David Radler, another Hollinger executive and a major shareholder in Ravelston, a total of $5.5 million in exchange for their promising not to compete with APC for three years after they stopped working for Hollinger. ….That Black and the others would start a newspaper in Mammoth Lake to compete with APC’s tiny newspaper there was ridiculous…. Although Hollinger is a large, sophisticated, public corporation, no document was found to indicate that the $5.5 million in payments was ever approved by the corporation or credited to the management-fees account on its books.

Black attempted to argue that the executives were entitled to the money as management fees, but the court made short shrift of these arguments, observing:

They are making a no harm-no foul argument, and such arguments usually fare badly in criminal cases. Suppose your employer owes you $100 but balks at paying, so you help yourself to the money from the cash register. That is theft, e.g., State v. Winston, 295 S.E.2d 46, 51 (W. Va. 1982); Edwards v. State, 181 N.W.2d 383, 387-88 (Wis. 1970); State v. Self, 713 P.2d 142, 144 (Wash. App. 1986), even though if the employer really owes you the money you have not harmed him. You are punishable because you are not entitled to take the law into your own hands. Harmlessness is rarely a defense to a criminal charge; if you embezzle money from your employer and replace it (with interest!) before the embezzlement is detected, you still are guilty of embezzlement.

The name of the newspaper caught my eye. Look at the following location map from Eli Rabett: Conrad Black’s Mammoth Lake CA (population 7093) is none other than the Mammoth Lakes CA (Wikipedia population 7093) of the location map, described at Wikipedia as follows: “Mammoth Lakes is located at 37°38′13″N, 118°58′44″W (37.636893, -118.978915), at an elevation of approximately 7900 feet (2400 m). As of the census[2] of 2000, there were 7,093 people, 2,814 households, and 1,516 families residing in the town.”

Maybe they should have characterized the fees as being for not collecting bristlecone pine data. The mind boggles at the possibilities.

The case attracted a lot of attention here because of Black’s using company property as though it were “personal” property. Data, and even correspondence, are also “property”, an issue that I’ll return to on another occasion.

UPDATE: Mammoth Lakes CA is extremely close to the sites studied by Constance Millar in Millar et al 2006 (Whitewing Mt, San Joaquin Ridge) in a study discussed on several occasions e.g. here, here. Millar says:

Preserved vegetation in the Mammoth Lakes – Long Valley Caldera region provides a detailed record of late Holocene climatic, vegetation, and volcanic dynamics. Scattered across the otherwise barren summits of Whitewing Mtn and San Joaquin Ridge, adjacent to and along the Sierra Nevada crest (Fig. 1), are abundant well-preserved deadwood tree stems (Fig. 2).

As noted on prior occasions, they date the destruction of the trees on Whitewing Mtn to “late summer AD1350” and model climate in the MWP as significantly warmer than at present (notwithstanding the Mannian PC1 which records near glacial coldness during this period in California).

Deadwood tree stems scattered above treeline on tephra-covered slopes of Whitewing Mtn (3051 m) and San Joaquin Ridge (3122 m) show evidence of being killed in an eruption from adjacent Glass Creek Vent, Inyo Craters. Using tree-ring methods, we dated deadwood to AD 815–1350 and infer from death dates that the eruption occurred in late summer AD 1350. Based on wood anatomy, we identified deadwood species as Pinus albicaulis, P. monticola, P. lambertiana, P. contorta, P. jeffreyi, and Tsuga mertensiana. Only P. albicaulis grows at these elevations currently; P. lambertiana is not locally native. Using contemporary distributions of the species, we modeled paleoclimate during the time of sympatry to be significantly warmer (+3.2°C annual minimum temperature) and slightly drier (−24 mm annual precipitation) than present, resembling values projected for California in the next 70–100 yr.

Here are location maps synchronizing the location map in Millar et al 2006 with a Google Earth image locating Mammoth Lakes CA on a somewhat similar scale.

 whitew29.gif  whitew30.jpg
By Stephen McIntyre | Posted in General | Tagged , | Comments (37)

Code: The Hockey Stick

Jul 1, 2008 – 2:55 PM


%%%% Matlab version of Hockey Stick
%
%Mann, M. E., R. S. Bradley, and M. K. Hughes (1998), Global-scale
%temperature patterns and climate forcing over the past six centuries,
%Nature, 392, 779– 787.
%
%Mann, M. E., R. S. Bradley, and M. K. Hughes (1999), Northern
%Hemisphere temperatures during the past millennium: Inferences, uncertainties,
%and limitations, Geophys. Res. Lett., 26, 759– 762.
%
% Tested with Matlab version 7.4. (R2007a)
% UC - uc_edit at yahoo.com - http://signals.auditblogs.com/
%%%%% hockeystick.txt (.m) ver 1.1 July 08

close all
%%% PART I %%%
%%% Get the data from web and write to files (once)%%%
%%% Your firewall might have a word at this point..%

if ~exist('downloaded_all.mat','file')

% Original reconstruction
urlwrite('http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/RECONS/nhem-recon.dat','nhem_recon.dat')

% MBH98 rec and Instrumental mean series

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt','nhmean.txt')

%
%urlwrite('ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/by_contributor/mann1998/nhem-dense.dat','nhem_dense.dat');

% Sparse instrumental data

urlwrite('http://www.ncdc.noaa.gov/paleo/ei/ei_data/nhem-sparse.dat','nhem_sparse.dat')

% gridpoints

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/gridpoints.txt','gridpoints.txt')

% proxies

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1400.txt','data1400.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1450.txt','data1450.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1500.txt','data1500.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1600.txt','data1600.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1700.txt','data1700.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1730.txt','data1730.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1750.txt','data1750.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1760.txt','data1760.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1780.txt','data1780.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1800.txt','data1800.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/PROXY/data1820.txt','data1820.txt')

% MBH99 AD1000

urlwrite('http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/PROXIES/pc1-fixed.dat','pc1_fixed.dat');

% http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/Millennium/DATA/PROXIES/
% morc014 is missing, so using CA version

urlwrite('http://www.climateaudit.info/data/mbh99/proxy.txt','proxy.txt')

% Eigenvalues for S-matrix
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/tpca-eigenvals.txt','tpca_eigenvals.txt')

% U matrix

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc01.txt','pc01.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc02.txt','pc02.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc03.txt','pc03.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc04.txt','pc04.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc05.txt','pc05.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc06.txt','pc06.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc07.txt','pc07.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc08.txt','pc08.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc09.txt','pc09.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc10.txt','pc10.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc11.txt','pc11.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc12.txt','pc12.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc13.txt','pc13.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc14.txt','pc14.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc15.txt','pc15.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/pc16.txt','pc16.txt')

% V matrix

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof01.txt','eof01.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof02.txt','eof02.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof03.txt','eof03.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof04.txt','eof04.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof05.txt','eof05.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof06.txt','eof06.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof07.txt','eof07.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof08.txt','eof08.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof09.txt','eof09.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof10.txt','eof10.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof11.txt','eof11.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof12.txt','eof12.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof13.txt','eof13.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof14.txt','eof14.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof15.txt','eof15.txt');
urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/eof16.txt','eof16.txt')

% Standard deviations

urlwrite('http://www.nature.com/nature/journal/v430/n6995/extref/Instrume/tpca-standard.txt','tpca_standard.txt')

downloaded_all=1;

save downloaded_all downloaded_all % just to mark successful download

end % if exist

%%% PART II %%%
%%% Load the data from files and preprocess for calibration%%%

load nhem_recon.dat

% Proxies

fid1=fopen('proxy.txt','r');
dds=[];
tline=fgetl(fid1);
while 1
tline = fgetl(fid1);
if ~ischar(tline), break, end

if tline(1:3)=='976', break,end % NAs

dds=[dds;str2num(tline)];

end

temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values (kind of manually..)
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:14)); temp2=str2num(tline(19:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

tline=fgetl(fid1);
temp=str2num(tline(1:13)); temp2=str2num(tline(18:end)); % replace NAs with previous values
dds=[dds; temp 1.063 temp2];

fclose(fid1);

ad1000=(dds(:,3:16));

ad1000b=ad1000; % fixed PC1 version
load pc1_fixed.dat
ad1000b(:,3)=pc1_fixed(:,2);

% standardize by detrended std ( http://www.nature.com/nature/journal/v430/n6995/extref/METHODS/README.txt )

AD1000=(ad1000-repmat(mean(ad1000(903:end,:)),981,1))./repmat(std(detrend(ad1000(903:end,:))),981,1);
AD1000b=(ad1000b-repmat(mean(ad1000b(903:end,:)),981,1))./repmat(std(detrend(ad1000b(903:end,:))),981,1);

load data1400.txt

ad1400=[zeros(400,22); data1400(:,2:end)]; % fix size to 981, helps later

for ii=1:981
iNaN=isnan(ad1400(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1400(ii,fiNaN)=ad1400(ii-1,fiNaN);
end
end

AD1400=(ad1400-repmat(mean(ad1400(903:end,:)),981,1))./repmat(std(detrend(ad1400(903:end,:))),981,1);

load data1450.txt
ad1450=[zeros(450,25); data1450(:,2:end)];

for ii=1:981
iNaN=isnan(ad1450(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1450(ii,fiNaN)=ad1450(ii-1,fiNaN);
end
end

AD1450=(ad1450-repmat(mean(ad1450(903:end,:)),981,1))./repmat(std(detrend(ad1450(903:end,:))),981,1);

load data1500.txt
ad1500=[zeros(500,28); data1500(:,2:end)];

for ii=1:981
iNaN=isnan(ad1500(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1500(ii,fiNaN)=ad1500(ii-1,fiNaN);
end
end

AD1500=(ad1500-repmat(mean(ad1500(903:end,:)),981,1))./repmat(std(detrend(ad1500(903:end,:))),981,1);

load data1600.txt
ad1600=[zeros(600,57); data1600(:,2:end)];

for ii=1:981
iNaN=isnan(ad1600(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1600(ii,fiNaN)=ad1600(ii-1,fiNaN);
end
end

AD1600=(ad1600-repmat(mean(ad1600(903:end,:)),981,1))./repmat(std(detrend(ad1600(903:end,:))),981,1);

load data1700.txt
ad1700=[zeros(700,74); data1700(:,2:end)];

for ii=1:981
iNaN=isnan(ad1700(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1700(ii,fiNaN)=ad1700(ii-1,fiNaN);
end
end

AD1700=(ad1700-repmat(mean(ad1700(903:end,:)),981,1))./repmat(std(detrend(ad1700(903:end,:))),981,1);

load data1730.txt
ad1730=[zeros(730,79); data1730(:,2:end)];

for ii=1:981
iNaN=isnan(ad1730(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1730(ii,fiNaN)=ad1730(ii-1,fiNaN);
end
end

AD1730=(ad1730-repmat(mean(ad1730(903:end,:)),981,1))./repmat(std(detrend(ad1730(903:end,:))),981,1);

load data1750.txt
ad1750=[zeros(750,89); data1750(:,2:end)];

for ii=1:981
iNaN=isnan(ad1750(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1750(ii,fiNaN)=ad1750(ii-1,fiNaN);
end
end

AD1750=(ad1750-repmat(mean(ad1750(903:end,:)),981,1))./repmat(std(detrend(ad1750(903:end,:))),981,1);

load data1760.txt
ad1760=[zeros(760,93); data1760(:,2:end)];
for ii=1:981
iNaN=isnan(ad1760(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1760(ii,fiNaN)=ad1760(ii-1,fiNaN);
end
end

AD1760=(ad1760-repmat(mean(ad1760(903:end,:)),981,1))./repmat(std(detrend(ad1760(903:end,:))),981,1);

load data1780.txt
ad1780=[zeros(780,97); data1780(:,2:end)];

for ii=1:981
iNaN=isnan(ad1780(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1780(ii,fiNaN)=ad1780(ii-1,fiNaN);
end
end

AD1780=(ad1780-repmat(mean(ad1780(903:end,:)),981,1))./repmat(std(detrend(ad1780(903:end,:))),981,1);

load data1800.txt
ad1800=[zeros(800,102); data1800(:,2:end)];
for ii=1:981
iNaN=isnan(ad1800(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1800(ii,fiNaN)=ad1800(ii-1,fiNaN);
end
end

AD1800=(ad1800-repmat(mean(ad1800(903:end,:)),981,1))./repmat(std(detrend(ad1800(903:end,:))),981,1);

load data1820.txt
ad1820=[zeros(820,112); data1820(:,2:end)];

for ii=1:981
iNaN=isnan(ad1820(ii,:)); %% if there are NaNs, replace with previous
if sum(iNaN)
fiNaN=find(iNaN);
ad1820(ii,fiNaN)=ad1820(ii-1,fiNaN);
end
end

AD1820=(ad1820-repmat(mean(ad1820(903:end,:)),981,1))./repmat(std(detrend(ad1820(903:end,:))),981,1);

% Instrumental data singular value decomposition, A=U*S*Vt

load pc01.txt
load pc02.txt
load pc03.txt
load pc04.txt
load pc05.txt
load pc06.txt
load pc07.txt
load pc08.txt
load pc09.txt
load pc10.txt
load pc11.txt
load pc12.txt
load pc13.txt
load pc14.txt
load pc15.txt
load pc16.txt

U=[pc01(:,2) pc02(:,2) pc03(:,2) pc04(:,2) pc05(:,2) pc06(:,2) pc07(:,2) pc08(:,2) pc09(:,2) pc10(:,2) pc11(:,2) pc12(:,2) pc13(:,2) pc14(:,2) pc15(:,2) pc16(:,2)];

load tpca_eigenvals.txt

S=diag(tpca_eigenvals(1:16,2)); % First 16 values

load eof01.txt
load eof02.txt
load eof03.txt
load eof04.txt
load eof05.txt
load eof06.txt
load eof07.txt
load eof08.txt
load eof09.txt
load eof10.txt
load eof11.txt
load eof12.txt
load eof13.txt
load eof14.txt
load eof15.txt
load eof16.txt

Vt=[eof01(:,2) eof02(:,2) eof03(:,2) eof04(:,2) eof05(:,2) eof06(:,2) eof07(:,2) eof08(:,2) eof09(:,2) eof10(:,2) eof11(:,2) eof12(:,2) eof13(:,2) eof14(:,2) eof15(:,2) eof16(:,2)]';

load gridpoints.txt

i_nh=find(ge(gridpoints(:,4),0)); % indices for NH in i_nh

% stds from tpca_standard.txt
load tpca_standard.txt
sigmas=tpca_standard(i_nh,3);

% Reference temperature

fid1=fopen('nhmean.txt','r');
dds=[];
tline=fgetl(fid1);
while 1
tline = fgetl(fid1);
if ~ischar(tline), break, end
dds=[dds;tline];
end

fclose(fid1);

ddn=str2num(dds);

RefT=ddn(503:581,3);

% Sparse reference

fid1=fopen('nhem_sparse.dat','r');
dds=[];
for iz=1:34
tline=fgetl(fid1);
end
while 1
tline = fgetl(fid1);
if ~ischar(tline), break, end
dds=[dds;str2num(tline)];
end

fclose(fid1);

RefT_ver=dds(1:48,2); % Verification reference 1854-1901, sparse grid

SRefT=dds(49:127,2); % Sparse for calibration period

i_ver=855:902; % Verification indices for reconstruction
i_cal=903:981; % Calibration indices

% cosines of latitudes (Cos) , i_nh indices to NH
Cos=cos(gridpoints(i_nh,4)*pi/180);

t=(1000:1980)'; % reconstruction years

% Reconstruction stats, see http://www.climateaudit.org/?p=647
RE=[]; % calibration REs
RE_ver=[];
OUT=[]; % output
s2CAL=[];
R2=[];
R2_ver=[];
CE=[];
s2VER=[];

%%% PART III %%%
%%% Perform the calibration in 12 steps %%%

% Plot the original reconstruction
figure(1)
subplot(211)
plot(t,nhem_recon(:,2))
hold on
title('Original in blue, PC1 not fixed ')
xlabel('Year')
ylabel('Temperature')

subplot(212)
plot(t,nhem_recon(:,2))
hold on
title('Original in blue, PC1 fixed ')
xlabel('Year')
ylabel('Temperature')

for step_i=1:13 %% Run the calibration step by step, last for the pc1_fixed

if step_i==1

Pmc=AD1000; % calibration detrended variance =1, centered to cal mean

Usel=[1]
chkd=1:400;

elseif step_i==2

chkd=401:450;
Pmc=AD1400;

elseif step_i==3

chkd=451:500;
Usel=[1 2] %% select U-vectors using http://holocene.meteo.psu.edu/shared/research/ONLINE-PREPRINTS/MultiProxy/stats-supp.html
Pmc=AD1450;

elseif step_i==4

chkd=501:600;
Usel=[1 2]
Pmc=AD1500;

elseif step_i==5

chkd=601:700;
Usel=[1 2 11 15]
Pmc=AD1600;

elseif step_i==6

chkd=701:730;
Usel=[1 2 5 11 15]
Pmc=AD1700;

elseif step_i==7

chkd=731:750;
Usel=[1 2 5 11 15]
Pmc=AD1730;

elseif step_i==8

chkd=751:760;
Usel=[1 2 3 5 6 8 11 15]
Pmc=AD1750;

elseif step_i==9

chkd=761:780;
Usel=[1 2 3 4 5 7 9 11 15]
Pmc=AD1760;

elseif step_i==10

chkd=781:800;
Usel=[1 2 3 4 5 7 9 11 14 15 16]
Pmc=AD1780;

elseif step_i==11

chkd=801:820;
Usel=[1 2 3 4 5 7 9 11 14 15 16]
Pmc=AD1800;

elseif step_i==12

chkd=821:981;
Usel=[1 2 3 4 5 7 9 11 14 15 16]

% max RE with 11 Us
%Usel=[ 1 2 3 4 5 7 8 9 11 12 13]

Pmc=AD1820;

else

Pmc=AD1000b; %
Usel=[1]
chkd=1:400;

end % if step_i

% http://www.climateaudit.org/?p=2445
zz=pinv(U(1:79,Usel))*Pmc(903:end,:); % part one of CCE, temperature PCs as targets
zzz=Pmc(1:end,:)*pinv(zz); % part two of CCE

Uhscaled=zeros(length(zzz),length(Usel));

for zii=1:length(Usel) % variance matching,

Uhscaled(:,zii)=zzz(:,zii)/std(zzz(903:end,zii))*std(U(1:79,Usel(zii)));

%Uhscaled(:,zii)=zzz(:,zii) % Don't touch the variance

end

temp_diag=diag(S);
temp_diag=temp_diag(Usel);
S_=diag(temp_diag); % Selected values of S diagonal

% http://www.climateaudit.org/?p=520
%
% http://www.climateaudit.org/?p=530

Recs=Uhscaled*S_*Vt(Usel,:); % back to temperatures

RecsN=Recs(:,i_nh)*diag(sigmas)*ones(length(i_nh),1)/sum(Cos); % ?, but works

RecNH=RecsN; %

V1=nhem_recon(chkd,2);
V2=RecNH(chkd,1);

figure(1)
if lt(step_i,13)

OUT=[OUT; t(chkd) V2 V2]; % OUT= [time non-fixed fixed]

subplot(211)
plot(t(chkd),V2,'g')

if gt(step_i,1)
subplot(212)
plot(t(chkd),V2,'g')
end
else
subplot(212)
plot(t(chkd),V2,'g')
OUT(1:400,3)=V2;
end

Residual=RecNH(903:981)-RefT;

Residual_ver=RecNH(i_ver)-RefT_ver;

s2CAL=[s2CAL; t(chkd(1)) 2*std(Residual)]; % Residual zero mean

s2VER=[s2VER; 2*sqrt(Residual_ver'*Residual_ver/length(Residual_ver)) ];

RE=[RE; t(chkd(1)) 1-Residual'*Residual/(RefT'*RefT)];

ctemp=corrcoef(RecNH(903:981),RefT);

R2=[R2; t(chkd(1)) ctemp(2,1)^2];

RE_ver=[RE_ver; t(chkd(1)) 1-Residual_ver'*Residual_ver/(RefT_ver'*RefT_ver)]; % mean(RefT)=0

CE=[CE; 1-Residual_ver'*Residual_ver/((RefT_ver-mean(RefT_ver))'*(RefT_ver-mean(RefT_ver))) ];

ctemp=corrcoef(RecNH(i_ver),RefT_ver);

R2_ver=[R2_ver; t(chkd(1)) ctemp(2,1)^2 ];

end % for i=

format long g

disp('Calibration R2: ')
R2

disp('Verification R2: ')
R2_ver

% plot the AD 1000-1850 trends

XX=[ones(851,1) (1000:1850)'];
Tnf=pinv(XX)*OUT(1:851,2);
Tf=pinv(XX)*OUT(1:851,3);

% trends per century
Tnf_c=Tnf(2)*100;
Tf_c=Tf(2)*100;

Trendnf=XX*Tnf;
Trendf=XX*Tf;

figure(1)
subplot(211)
plot((1000:1850)',Trendnf,'r')
text(1500,0.25,['Trend per century: ', num2str(Tnf_c,3) ])
subplot(212)
plot((1000:1850)',Trendf,'r')
text(1500,0.25,['Trend per century: ', num2str(Tf_c,3) ])

% That's all


Steve: also R-analysis http://www.climateaudit.org/tag/algebra

By uc00 | Posted in MBH98, Replication, Source Code, Uncategorized | Tagged , , , | Comments Off on Code: The Hockey Stick

Sea Ice – May 2008 Results

Jun 30, 2008 – 7:25 PM

Sea ice continues to get lots of attention. I previously discussed March 2008 results here. June 2008 results should be available soon, but in the mean time, I’ve updated my own graphics showing a tripartite image of global, SH and NH, instead of only showing NH sea ice, as done in other recent comments e.g. Phil here, the UK Independent and Michael Tobis here. Continue reading →

By Stephen McIntyre | Posted in General | Comments (132)