clear
echo on
%Solution to Problem 11
%
%This assignment asked you to find the
%average global temperature over the last
%5 centuries off of the net, and to do a
%regression analysis on the data from
%1900 to the present.  The data could
%be found using the altavista search
%engine (http://altavista.digital.com)
%and doing an advanced keyword search.
%The search I used was:
%
%global near temperature and 1434 and 1922 and 1933
%
%which returned one record.
%
%The URL we wanted was:
%
%www.iphc.washington.edu/staff/hare/html/decadal/post1977/nhtemp.txt
%
%which is from a paper in Nature.  Search of the 
%web using the parameters:
%
%Mann near Bradley near Hughes
%
%for the authors of the paper yields a large number of matches - it
%was a very popular article!  On the first page of results, however
%is a link to Mann's home page at the University of Virginia. A
%couple of clicks later you get to the url:
%
%ftp://eclogite.geo.umass.edu/pub/mann/ONLINE-PREPRINTS/MultiProxy/mbh98.pdf
%
%which provides a .pdf (Adobe portable document format) file containing
%the entire article, which you can print out.
%
%If you print the article out, you find
%that the data table isn't in the article itself, although
%the plot that you get from the data is the same as that
%in the paper.  It is not clear the relationship between
%the original website and the authors of the paper.  The website
%is actually on the server of the International Pacific
%Halibut Commission! 
%
%It is always important to know the source of the
%data you get off of the internet.  There isn't much
%peer review in cyberspace!
pause
%The data is attached below this document.  Because the
%data is not in a simple matrix form, matlab can't just
%read it in as one piece.  To deal with this we will break
%it into three parts, saving the first part under the name
%early.dat, the second under the name middle.dat, and the
%third under the name late.dat.  We can then read each of
%these files in, and concatenate the appropriate columns to
%get the complete time series.  This saves -alot- of typing!
pause
%Now we load the data in (we presume you have saved the files
%already!):
load 'early.dat'
load 'middle.dat'
load 'late.dat'
%and we operate to get the two time series of interest:
timerecon=[early(:,1);middle(:,1)];
temprecon=[early(:,2);middle(:,3)];
timeactual=[middle(:,1);late(:,1)];
tempactual=[middle(:,2);late(:,2)];
%and we plot the data up:
figure(1)
plot(timerecon,temprecon,'og',timeactual,tempactual,'xr')
xlabel('year')
ylabel('temperature anomaly (deg. C)')
labels=str2mat('reconstruction','measured');
legend(labels)
pause
%Note that the plot is of the temperature anomaly -
%the deviation from the baseline temperature which
%(if you read the paper) is the average over the range
%1902-1980 (the middle time period).
pause
%OK, now for the regression part.  The part of interest
%is the measured temperature anomaly.  We want to fit
%this data to a line.  Thus:
n=length(timeactual);
a=[ones(n,1),timeactual];
k=inv(a'*a)*a';
x=k*tempactual
figure(1)
hold on
plot(timeactual,a*x,'r')
hold off
labels=str2mat(labels,'regression fit')
legend(labels)
pause
%Note that the average temperature rise per year is really
%very small!  This coefficient is only:
x(2)
%which is imperceptible - far less than the day-to-day variation.
pause
%We are also asked to determine how much the temperature
%should go up between now and the year 2020.  This is just:
tempinc=x(2)*20
%which is also pretty small!  This prediction can be added to
%the plot as well:
figure(1)
hold on
timefit=[1995;2020];
afit=[ones(2,1),timefit];
plot(timefit,afit*x,'b')
plot(2020,[1,2020]*x,'c*')
labels=str2mat(labels,'extrapolation','value at 2020');
legend(labels)
hold off
pause
%We can also calculate the error in the
%regression parameters by looking at the
%scatter in the data.  We will be discussing
%this in class in the next few lectures.  We
%begin with the residual:

format short e
r=a*x-tempactual;
vartemp=r'*r/(n-2)
%which is the population variance of the
%individual annual temperature measure-
%ments.
pause

%We also have the matrix of covariance of the
%fitting parameters:
varx=k*k'*vartemp

pause
%and finally, we have the rate of global
%warming given by x(2):
[x(2),varx(2,2)^.5]
%Where the error assumes that there is
%no covariance in the temperature measurements.
%Note that this is almost certainly -not-
%true, as cycles in the sun and other effects
%may persist for decades or centuries!  Still
%the regression analysis does suggest why
%global warming is considered significant!
echo off



________________

Here is the data.  Break it up and save it in three parts
under the appropriate names!

These data were published in the followng paper:

Mann, M.E., Bradley, R.S., Hughes, M.K. 1998.  Global-Scale Temperature Patterns 
and Climate Forcing Over the Past Six Centuries, Nature, 392, 779-787.

Here is the abstract to the paper:

Reconstructions of annual global surface temperature patterns over several
centuries are now possible, based on the multivariate calibration of widely 
distributed high-resolution proxy climate indicators.  These reconstructions 
provide insight into both the spatial and temporal nature of climatic 
variations during the past six centuries.  Time-dependent correlations of 
these temperature reconstructions with time series representing greenhouse 
gases, solar irradiance, and volcanic aerosols, suggest that each of these 
forcings has played a role in the climatic variability of the past several 
centuries, with greenhouse gases appearing to emerge as the dominant 
forcing during the 20th century. Northern hemisphere mean annual temperatures 
for three of the past eight years are warmer than any other year since 
(at least) 1400 AD, at a greater than 99.5\% level of confidence.


Year    Recon.  Instrum.
1400    -0.305  
1401    -0.239  
1402    -0.132  
1403    -0.092  
1404    -0.217  
1405    -0.152  
1406    -0.050  
1407    -0.113  
1408    -0.102  
1409    -0.086  
1410    -0.109  
1411    -0.076  
1412    -0.179  
1413    -0.157  
1414    -0.138  
1415    -0.090  
1416    -0.003  
1417    -0.081  
1418    -0.106  
1419    -0.147  
1420    -0.164  
1421    -0.286  
1422    -0.141  
1423    -0.080  
1424    -0.016  
1425    -0.127  
1426    -0.074  
1427    -0.208  
1428    -0.154  
1429    -0.182  
1430    -0.122  
1431    -0.201  
1432    -0.177  
1433    -0.046  
1434    -0.097  
1435    -0.032  
1436    -0.029  
1437    -0.009  
1438    -0.145  
1439    -0.176  
1440    -0.040  
1441    0.015   
1442    -0.093  
1443    -0.020  
1444    -0.030  
1445    -0.191  
1446    -0.194  
1447    -0.111  
1448    -0.188  
1449    -0.234  
1450    -0.251  
1451    -0.200  
1452    -0.149  
1453    -0.238  
1454    -0.317  
1455    -0.316  
1456    -0.354  
1457    -0.327  
1458    -0.413  
1459    -0.387  
1460    -0.342  
1461    -0.390  
1462    -0.459  
1463    -0.331  
1464    -0.456  
1465    -0.382  
1466    -0.339  
1467    -0.347  
1468    -0.400  
1469    -0.315  
1470    -0.282  
1471    -0.415  
1472    -0.400  
1473    -0.353  
1474    -0.153  
1475    -0.324  
1476    -0.241  
1477    -0.193  
1478    -0.200  
1479    -0.136  
1480    -0.152  
1481    -0.132  
1482    -0.256  
1483    -0.186  
1484    -0.091  
1485    -0.091  
1486    -0.173  
1487    -0.171  
1488    -0.106  
1489    -0.079  
1490    -0.100  
1491    -0.073  
1492    -0.272  
1493    -0.081  
1494    -0.130  
1495    -0.224  
1496    -0.308  
1497    -0.346  
1498    -0.233  
1499    -0.219  
1500    -0.191  
1501    -0.204  
1502    -0.259  
1503    -0.108  
1504    -0.096  
1505    -0.177  
1506    -0.173  
1507    -0.066  
1508    -0.043  
1509    -0.012  
1510    0.063   
1511    0.031   
1512    -0.056  
1513    -0.167  
1514    -0.123  
1515    -0.159  
1516    -0.074  
1517    -0.133  
1518    -0.263  
1519    -0.206  
1520    -0.144  
1521    -0.099  
1522    -0.232  
1523    -0.121  
1524    -0.107  
1525    -0.046  
1526    -0.088  
1527    -0.033  
1528    -0.195  
1529    -0.173  
1530    -0.045  
1531    -0.017  
1532    -0.077  
1533    -0.158  
1534    -0.161  
1535    -0.059  
1536    -0.112  
1537    -0.157  
1538    -0.196  
1539    -0.059  
1540    -0.144  
1541    -0.179  
1542    -0.112  
1543    -0.195  
1544    -0.172  
1545    -0.133  
1546    -0.067  
1547    -0.219  
1548    -0.210  
1549    -0.224  
1550    -0.165  
1551    -0.147  
1552    -0.162  
1553    -0.179  
1554    -0.052  
1555    -0.055  
1556    -0.030  
1557    -0.118  
1558    -0.199  
1559    -0.105  
1560    -0.136  
1561    -0.168  
1562    -0.091  
1563    0.004   
1564    0.039   
1565    -0.017  
1566    0.022   
1567    -0.109  
1568    -0.029  
1569    -0.067  
1570    -0.067  
1571    -0.081  
1572    -0.159  
1573    -0.173  
1574    -0.020  
1575    -0.126  
1576    -0.114  
1577    -0.224  
1578    -0.294  
1579    -0.295  
1580    -0.233  
1581    -0.176  
1582    -0.092  
1583    -0.129  
1584    -0.111  
1585    -0.126  
1586    -0.199  
1587    -0.282  
1588    -0.227  
1589    -0.239  
1590    -0.297  
1591    -0.298  
1592    -0.203  
1593    -0.160  
1594    -0.178  
1595    -0.224  
1596    -0.167  
1597    -0.148  
1598    -0.155  
1599    -0.155  
1600    -0.129  
1601    -0.366  
1602    -0.194  
1603    -0.250  
1604    -0.268  
1605    -0.168  
1606    -0.263  
1607    -0.199  
1608    -0.223  
1609    -0.152  
1610    -0.102  
1611    -0.232  
1612    -0.192  
1613    -0.124  
1614    -0.226  
1615    -0.208  
1616    -0.143  
1617    -0.063  
1618    -0.080  
1619    -0.035  
1620    -0.033  
1621    -0.075  
1622    -0.165  
1623    -0.249  
1624    -0.189  
1625    -0.248  
1626    -0.344  
1627    -0.252  
1628    -0.166  
1629    -0.179  
1630    -0.118  
1631    -0.260  
1632    -0.073  
1633    -0.156  
1634    -0.024  
1635    0.003   
1636    0.086   
1637    -0.089  
1638    -0.149  
1639    -0.138  
1640    -0.022  
1641    -0.281  
1642    -0.382  
1643    -0.193  
1644    -0.256  
1645    -0.207  
1646    -0.200  
1647    -0.206  
1648    -0.149  
1649    -0.164  
1650    -0.107  
1651    -0.118  
1652    -0.061  
1653    -0.082  
1654    -0.130  
1655    -0.040  
1656    -0.017  
1657    -0.024  
1658    -0.058  
1659    -0.100  
1660    -0.121  
1661    -0.095  
1662    -0.100  
1663    -0.196  
1664    -0.297  
1665    -0.258  
1666    -0.132  
1667    -0.329  
1668    -0.263  
1669    -0.305  
1670    -0.208  
1671    -0.167  
1672    -0.294  
1673    -0.383  
1674    -0.305  
1675    -0.293  
1676    -0.345  
1677    -0.195  
1678    -0.299  
1679    -0.298  
1680    -0.339  
1681    -0.125  
1682    -0.176  
1683    -0.264  
1684    -0.255  
1685    -0.278  
1686    -0.247  
1687    -0.164  
1688    -0.244  
1689    -0.072  
1690    -0.057  
1691    -0.072  
1692    -0.160  
1693    -0.095  
1694    -0.063  
1695    -0.225  
1696    -0.372  
1697    -0.259  
1698    -0.347  
1699    -0.362  
1700    -0.393  
1701    -0.315  
1702    -0.172  
1703    -0.189  
1704    -0.236  
1705    -0.362  
1706    -0.269  
1707    -0.145  
1708    -0.172  
1709    -0.196  
1710    -0.176  
1711    -0.237  
1712    -0.255  
1713    -0.211  
1714    -0.264  
1715    -0.187  
1716    -0.257  
1717    -0.241  
1718    -0.145  
1719    -0.161  
1720    -0.084  
1721    -0.096  
1722    -0.077  
1723    0.057   
1724    -0.028  
1725    -0.273  
1726    -0.063  
1727    -0.014  
1728    -0.190  
1729    -0.050  
1730    -0.102  
1731    -0.165  
1732    -0.252  
1733    -0.176  
1734    -0.257  
1735    -0.121  
1736    -0.113  
1737    -0.117  
1738    -0.141  
1739    -0.093  
1740    -0.082  
1741    -0.081  
1742    -0.225  
1743    -0.251  
1744    -0.275  
1745    -0.199  
1746    -0.255  
1747    -0.181  
1748    -0.058  
1749    -0.006  
1750    -0.121  
1751    -0.138  
1752    -0.145  
1753    -0.149  
1754    -0.116  
1755    -0.156  
1756    -0.032  
1757    -0.171  
1758    -0.140  
1759    -0.021  
1760    0.037   
1761    -0.144  
1762    -0.037  
1763    -0.166  
1764    -0.231  
1765    -0.076  
1766    -0.048  
1767    0.022   
1768    -0.063  
1769    -0.219  
1770    0.091   
1771    0.012   
1772    0.122   
1773    -0.041  
1774    0.044   
1775    0.052   
1776    -0.011  
1777    -0.112  
1778    -0.086  
1779    -0.184  
1780    -0.206  
1781    -0.059  
1782    -0.251  
1783    -0.091  
1784    -0.128  
1785    -0.267  
1786    -0.266  
1787    -0.262  
1788    -0.120  
1789    -0.166  
1790    -0.265  
1791    -0.062  
1792    -0.018  
1793    -0.078  
1794    0.027   
1795    -0.117  
1796    -0.098  
1797    -0.192  
1798    -0.097  
1799    -0.180  
1800    -0.076  
1801    -0.098  
1802    -0.168  
1803    -0.254  
1804    -0.032  
1805    -0.182  
1806    -0.050  
1807    -0.142  
1808    -0.079  
1809    -0.162  
1810    -0.131  
1811    -0.127  
1812    -0.239  
1813    -0.257  
1814    -0.274  
1815    -0.262  
1816    -0.325  
1817    -0.340  
1818    -0.259  
1819    -0.318  
1820    -0.442  
1821    -0.211  
1822    -0.217  
1823    -0.247  
1824    -0.269  
1825    -0.205  
1826    -0.055  
1827    -0.149  
1828    -0.109  
1829    -0.164  
1830    -0.140  
1831    -0.119  
1832    -0.237  
1833    -0.119  
1834    0.064   
1835    -0.307  
1836    -0.214  
1837    -0.370  
1838    -0.425  
1839    -0.216  
1840    -0.375  
1841    -0.198  
1842    -0.364  
1843    -0.245  
1844    -0.245  
1845    -0.312  
1846    -0.088  
1847    -0.276  
1848    -0.204  
1849    -0.297  
1850    -0.245  
1851    -0.287  
1852    -0.136  
1853    -0.227  
1854    -0.255  
1855    -0.170  
1856    -0.097  
1857    -0.160  
1858    -0.257  
1859    -0.126  
1860    -0.286  
1861    -0.195  
1862    -0.192  
1863    -0.137  
1864    -0.389  
1865    -0.140  
1866    -0.197  
1867    -0.203  
1868    -0.167  
1869    -0.233  
1870    -0.390  
1871    -0.328  
1872    -0.236  
1873    -0.209  
1874    -0.257  
1875    -0.236  
1876    -0.189  
1877    -0.014  
1878    -0.087  
1879    -0.296  
1880    -0.237  
1881    -0.198  
1882    -0.204  
1883    -0.249  
1884    -0.212  
1885    -0.190  
1886    -0.108  
1887    -0.326  
1888    -0.169  
1889    -0.134  
1890    -0.311  
1891    -0.175  
1892    -0.319  
1893    -0.324  
1894    -0.197  
1895    -0.158  
1896    -0.080  
1897    -0.054  
1898    -0.220  
1899    -0.349  
1900    -0.125  
1901    -0.158  
1902    -0.366  -0.204
1903    -0.235  -0.280
1904    -0.385  -0.351
1905    -0.199  -0.171
1906    -0.130  -0.104
1907    -0.346  -0.371
1908    -0.362  -0.360
1909    -0.411  -0.290
1910    -0.277  -0.292
1911    -0.212  -0.339
1912    -0.326  -0.344
1913    -0.229  -0.368
1914    -0.175  -0.121
1915    -0.107  0.062
1916    -0.196  -0.235
1917    -0.319  -0.392
1918    -0.213  -0.226
1919    -0.083  -0.269
1920    -0.187  -0.203
1921    -0.156  -0.112
1922    -0.280  -0.200
1923    -0.044  -0.168
1924    -0.068  -0.197
1925    -0.164  -0.026
1926    0.119   0.110
1927    -0.032  -0.004
1928    0.068   -0.014
1929    -0.170  -0.206
1930    0.153   0.090
1931    -0.017  0.134
1932    0.024   0.076
1933    -0.114  -0.121
1934    0.074   0.061
1935    0.133   -0.011
1936    0.107   0.130
1937    0.240   0.174
1938    0.188   0.230
1939    0.167   0.147
1940    0.003   0.067
1941    0.162   0.128
1942    0.104   0.118
1943    0.162   0.145
1944    0.278   0.293
1945    0.136   0.092
1946    0.136   0.155
1947    0.225   0.074
1948    0.165   0.153
1949    -0.004  0.125
1950    0.057   -0.010
1951    -0.020  0.167
1952    0.091   0.198
1953    0.028   0.313
1954    0.075   0.050
1955    0.114   0.101
1956    -0.066  -0.139
1957    0.129   0.199
1958    0.264   0.305
1959    0.189   0.203
1960    0.260   0.170
1961    0.165   0.286
1962    0.082   0.312
1963    0.138   0.293
1964    0.108   -0.009
1965    0.040   -0.026
1966    0.241   0.128
1967    0.102   0.145
1968    0.173   0.043
1969    0.189   0.053
1970    0.133   0.099
1971    -0.108  -0.058
1972    0.039   -0.063
1973    0.106   0.189
1974    -0.101  -0.089
1975    0.042   0.020
1976    0.046   -0.199
1977    0.169   0.189
1978    0.098   0.107
1979    0.257   0.226
1980    0.121   0.213
1981            0.382
1982            0.156
1983            0.399
1984            0.102
1985            0.039
1986            0.198
1987            0.359
1988            0.413
1989            0.379
1990            0.574
1991            0.424
1992            0.270
1993            0.300
1994            0.450
1995            0.660