%linear CO2-temperature anomaly analysis to get HIW %A.Jarvis %initiated - Jan. 2023 (following LEC380 sessions from Jan 2015 onwards) %last edited - Nov. 2024 %this edition is still annual based on the 2023 data %but adds in the monthly updated 2024 average %offered in good faith but needs a good tidy up! clear all %Law Dome ice core pCO2 data %https://cdiac.ess-dive.lbl.gov/ftp/trends/co2/lawdome.combined.dat %https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/law/law2006.xls %YEAR pCO2(ppm) (uncertainty quoted at +/- 1ppm undefined) M = [2006 378.7 2005 376.7 2004 374.7 2003 372.8 2002 370.5 2001 368.3 2000 366.8 1999 365.5 1998 363.6 1997 360.3 1996 359.8 1995 358.3 1994 356.6 1993 354.9 1992 353.8 1991 352.6 1990 350.3 1989 349.5 1988 348.3 1987 346.1 1986 344.5 1985 343.2 1984 341.5 1983 340.1 1982 338.1 1981 337.6 1980 336.6 1979 334.2 1978 333.5 1977 331.7 1976 331.2 1974 328.1 1973 327.8 1972 324.1 1971 325.0 1970 324.0 1969 323.7 1967 322.9 1966 319.1 1964 318.6 1963 318.0 1961 315.7 1959 316.3 1958 314.4 1957 314.0 1956 316.3 1955 314.0 1954 311.8 1953 312.1 1950 312.6 1949 310.8 1948 310.5 1947 310.7 1946 311.5 1945 309.7 1944 311.6 1943 310.8 1942 311.3 1941 310.6 1940 311.4 1939 310.7 1938 309.6 1937 308.4 1935 308.5 1933 307.5 1929 306.7 1928 305.2 1927 305.0 1925 304.4 1924 305.2 1923 303.2 1919 303.4 1916 301.3 1914 300.4 1913 300.7 1912 298.4 1909 300.4 1906 297.7 1905 299.0 1904 295.1 1902 295.3 1901 296.5 1899 295.6 1896 298.2 1894 293.8 1893 294.6 1892 294.7 1889 291.9 1887 293.7 1886 290.6 1884 289.4 1883 291.9 1878 288.8 1874 290.5 1873 287.2 1870 287.4 1869 287.7 1867 285.2 1864 285.4 1862 286.6 1859 286.5 1855 284.9 1854 287.0 1852 288.6 1851 285.2 1849 287.7 1848 286.1 1846 284.1 1844 286.5 1841 283.0 1838 284.1 1835 283.7 1833 284.5 1827 285.1 1826 281.3 1814 284.3 1800 281.1 1799 283.7 1796 281.6 1794 281.5 1781 276.8 1780 279.5 1774 277.8 1763 276.5 1752 276.8 1749 276.9 1743 276.7 1734 278.2 1723 277.2 1694 276.5 1690 276.3 1682 275.9 1649 277.2 1640 276.6 1629 274.5 1610 271.8 1603 274.3 1591 278.7 1588 281.0 1573 281.9 1560 281.7 1550 282.8 1530 283.2 1502 282.4 1469 279.6 1449 281.7 1431 282.5 1429 279.5 1411 279.6 1391 280.0 1390 280.4 1350 280.1 1330 283.4 1306 281.5 1276 281.1 1258 282.1 1246 281.7 1207 283.6 1193 283.9 1160 283.9 1137 283.8 1105 282.8 1088 282.4 1058 282.8 1037 280.3 1025 280.8 1005 279.9 968 278.5 944 279.1 897 278.9 857 279.3 799 278.5 764 278.5 730 278.5 698 279.7 668 279.4 632 278.3 596 276.9 572 277.6 537 276.0 500 276.4 461 276.7 428 276.9 365 277.0 329 278.9 302 279.8 274 280.1 228 281.5 202 280.7 168 280.1 136 278.1 104 277.5 56 277.4 30 277.9 13 276.7]; yearLawDome = flipud(round(M(:,1))); LawDomeCO2 = flipud(M(:,2)); %Annual Mauna Loa atmospheric pCO2 data to 2023 %https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_annmean_mlo.txt %YEAR [pCO2](ppm) 'UNCERTAINTY'(UNDEFINED?) M = [1959 315.98 0.12 1960 316.91 0.12 1961 317.64 0.12 1962 318.45 0.12 1963 318.99 0.12 1964 319.62 0.12 1965 320.04 0.12 1966 321.37 0.12 1967 322.18 0.12 1968 323.05 0.12 1969 324.62 0.12 1970 325.68 0.12 1971 326.32 0.12 1972 327.46 0.12 1973 329.68 0.12 1974 330.19 0.12 1975 331.13 0.12 1976 332.03 0.12 1977 333.84 0.12 1978 335.41 0.12 1979 336.84 0.12 1980 338.76 0.12 1981 340.12 0.12 1982 341.48 0.12 1983 343.15 0.12 1984 344.87 0.12 1985 346.35 0.12 1986 347.61 0.12 1987 349.31 0.12 1988 351.69 0.12 1989 353.20 0.12 1990 354.45 0.12 1991 355.70 0.12 1992 356.54 0.12 1993 357.21 0.12 1994 358.96 0.12 1995 360.97 0.12 1996 362.74 0.12 1997 363.88 0.12 1998 366.84 0.12 1999 368.54 0.12 2000 369.71 0.12 2001 371.32 0.12 2002 373.45 0.12 2003 375.98 0.12 2004 377.70 0.12 2005 379.98 0.12 2006 382.09 0.12 2007 384.02 0.12 2008 385.83 0.12 2009 387.64 0.12 2010 390.10 0.12 2011 391.85 0.12 2012 394.06 0.12 2013 396.74 0.12 2014 398.81 0.12 2015 401.01 0.12 2016 404.41 0.12 2017 406.76 0.12 2018 408.72 0.12 2019 411.65 0.12 2020 414.21 0.12 2021 416.41 0.12 2022 418.53 0.12 2023 421.08 0.12]; yearMaunaLoa = M(:,1); MaunaLoaCO2 = M(:,2); %add 2024 estimate %using the NOAA de-seasoned trend estimate %and one month lagged to match the HadCRUT5 release frequency %https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt M = [2024 1 2024.0417 422.80 422.52 27 0.70 0.26 2024 2 2024.1250 424.55 423.63 22 1.24 0.51 2024 3 2024.2083 425.38 423.95 22 0.99 0.40 2024 4 2024.2917 426.51 424.01 22 0.99 0.40 2024 5 2024.3750 426.90 423.64 29 0.76 0.27 2024 6 2024.4583 426.91 424.47 20 0.65 0.28 2024 7 2024.5417 425.55 425.10 24 0.69 0.27 2024 8 2024.6250 422.99 424.82 22 1.08 0.44 2024 9 2024.7083 422.03 425.42 18 0.41 0.18 2024 10 2024.7917 422.38 425.65 22 0.35 0.14]; yearMaunaLoa = [yearMaunaLoa;2024]; MaunaLoaCO2 = [MaunaLoaCO2;M(end-1,5)]; %HadCRUT5A global temperature anomaly to 2023 %https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/analysis/diagnostics/HadCRUT.5.0.1.0.analysis.summary_series.global.annual.csv %YEAR HadCRUT5A(degC) 5thPERCENTILE(degC) 95thPERCENTILE(degC) M = [1850 -0.4177114 -0.58925647 -0.2461663 1851 -0.2333498 -0.41186792 -0.054831687 1852 -0.22939907 -0.40938243 -0.04941572 1853 -0.27035445 -0.43000934 -0.110699534 1854 -0.29152083 -0.4327115 -0.15033019 1855 -0.29691675 -0.43931878 -0.15451474 1856 -0.32035372 -0.46809322 -0.1726142 1857 -0.46723005 -0.61632216 -0.31813794 1858 -0.3887657 -0.53688604 -0.24064532 1859 -0.28126517 -0.42392084 -0.1386095 1860 -0.39016518 -0.5389766 -0.24135375 1861 -0.42911294 -0.597001 -0.26122484 1862 -0.5363694 -0.70365995 -0.3690788 1863 -0.34424406 -0.53403676 -0.15445141 1864 -0.46546507 -0.64809304 -0.28283712 1865 -0.33248132 -0.5244569 -0.14050578 1866 -0.3412875 -0.5218009 -0.16077414 1867 -0.35699412 -0.5530232 -0.16096504 1868 -0.35182714 -0.5294787 -0.1741756 1869 -0.31659195 -0.47646335 -0.15672056 1870 -0.32792753 -0.46869946 -0.18715557 1871 -0.36856276 -0.51413316 -0.22299235 1872 -0.32811058 -0.46317002 -0.19305114 1873 -0.3412969 -0.47250164 -0.21009217 1874 -0.3732512 -0.5071426 -0.2393598 1875 -0.37562594 -0.514041 -0.23721085 1876 -0.4241099 -0.56287116 -0.28534868 1877 -0.10110884 -0.22982001 0.02760234 1878 -0.011315192 -0.13121258 0.1085822 1879 -0.30363432 -0.43406436 -0.1732043 1880 -0.31583205 -0.44015092 -0.19151321 1881 -0.23224552 -0.35793498 -0.10655606 1882 -0.29553008 -0.4201501 -0.17091006 1883 -0.3464744 -0.4608177 -0.23213112 1884 -0.49232006 -0.6026686 -0.38197154 1885 -0.47112358 -0.5830682 -0.35917896 1886 -0.42090362 -0.5225382 -0.31926903 1887 -0.49878576 -0.61655986 -0.38101166 1888 -0.37937889 -0.49332377 -0.265434 1889 -0.24989556 -0.37222093 -0.12757017 1890 -0.50685817 -0.6324095 -0.3813068 1891 -0.40131494 -0.5373699 -0.26525998 1892 -0.5075585 -0.64432853 -0.3707885 1893 -0.49461922 -0.6315314 -0.35770702 1894 -0.48376393 -0.6255681 -0.34195974 1895 -0.4487516 -0.58202064 -0.3154826 1896 -0.28400728 -0.41740146 -0.15061307 1897 -0.25980017 -0.39852425 -0.12107606 1898 -0.48579213 -0.6176492 -0.35393503 1899 -0.35543364 -0.48639694 -0.22447036 1900 -0.2344939 -0.3664699 -0.1025179 1901 -0.29341024 -0.42826617 -0.15855433 1902 -0.43895653 -0.57330817 -0.3046049 1903 -0.5332871 -0.65805364 -0.40852055 1904 -0.59751105 -0.72495615 -0.47006598 1905 -0.40779322 -0.53042954 -0.2851569 1906 -0.31910878 -0.4449256 -0.19329198 1907 -0.5040763 -0.61996025 -0.38819233 1908 -0.5138197 -0.6301149 -0.39752448 1909 -0.53568715 -0.64385575 -0.42751858 1910 -0.5309095 -0.6454714 -0.4163476 1911 -0.539079 -0.6510565 -0.42710155 1912 -0.47553864 -0.5783255 -0.37275177 1913 -0.4670111 -0.5772189 -0.3568033 1914 -0.26243657 -0.37001243 -0.15486072 1915 -0.19167219 -0.30920455 -0.07413985 1916 -0.42002314 -0.5453339 -0.2947124 1917 -0.5428197 -0.67724293 -0.40839645 1918 -0.4243641 -0.56627965 -0.28244856 1919 -0.32528907 -0.46366385 -0.18691428 1920 -0.29835507 -0.4310807 -0.16562946 1921 -0.24044435 -0.36450884 -0.11637985 1922 -0.3390137 -0.449683 -0.22834438 1923 -0.31768188 -0.4260641 -0.20929964 1924 -0.3118017 -0.4204113 -0.2031921 1925 -0.28214198 -0.39495143 -0.16933253 1926 -0.122555 -0.23084709 -0.014262918 1927 -0.2291136 -0.33219776 -0.12602946 1928 -0.20646581 -0.31743687 -0.09549474 1929 -0.39244303 -0.50286 -0.28202602 1930 -0.17680542 -0.29041144 -0.06319938 1931 -0.103397675 -0.2126916 0.005896252 1932 -0.14546171 -0.25195518 -0.03896824 1933 -0.32234442 -0.4271004 -0.21758842 1934 -0.17433685 -0.27400395 -0.07466974 1935 -0.20605923 -0.30349734 -0.10862113 1936 -0.16952094 -0.26351926 -0.07552263 1937 -0.019198947 -0.11975877 0.08136088 1938 -0.012200737 -0.110303745 0.08590227 1939 -0.040797204 -0.14670469 0.06511029 1940 0.07593582 -0.041949686 0.19382133 1941 0.0381293 -0.16225392 0.23851252 1942 0.001406068 -0.19521242 0.19802456 1943 0.006421582 -0.19958311 0.21242628 1944 0.14410514 -0.05449483 0.3427051 1945 0.043088354 -0.1572829 0.24345961 1946 -0.11881461 -0.26595104 0.028321814 1947 -0.09120561 -0.23179047 0.049379256 1948 -0.12466127 -0.25913337 0.009810847 1949 -0.14380223 -0.25407746 -0.033526964 1950 -0.2266218 -0.33265698 -0.120586626 1951 -0.06115396 -0.15035023 0.028042309 1952 0.015354548 -0.08293598 0.11364508 1953 0.07763075 -0.020529611 0.17579111 1954 -0.11675023 -0.20850272 -0.02499773 1955 -0.19730994 -0.28442997 -0.11018991 1956 -0.26316562 -0.33912566 -0.18720558 1957 -0.035334915 -0.10056861 0.029898778 1958 -0.017632563 -0.08307456 0.04780944 1959 -0.048004813 -0.11036374 0.014354111 1960 -0.11545958 -0.17413756 -0.056781605 1961 -0.019999769 -0.07077983 0.030780297 1962 -0.06404272 -0.117309235 -0.010776198 1963 -0.036810614 -0.09057977 0.016958548 1964 -0.30586153 -0.34949452 -0.26222858 1965 -0.20442048 -0.2536107 -0.15523025 1966 -0.1488976 -0.19840918 -0.09938603 1967 -0.117539294 -0.1606568 -0.07442179 1968 -0.16864756 -0.21327549 -0.12401963 1969 -0.03138624 -0.071891375 0.009118891 1970 -0.08506408 -0.12605265 -0.044075526 1971 -0.20588905 -0.24447224 -0.16730586 1972 -0.09379131 -0.1316835 -0.055899136 1973 0.04995016 0.013479018 0.0864213 1974 -0.17252657 -0.21022564 -0.13482748 1975 -0.110754214 -0.15130511 -0.07020333 1976 -0.2158369 -0.25586063 -0.17581315 1977 0.1030885 0.06005668 0.14612031 1978 0.005255972 -0.034596615 0.045108557 1979 0.09085814 0.062358625 0.11935765 1980 0.19607204 0.16280398 0.2293401 1981 0.25001204 0.21939126 0.28063282 1982 0.034268282 -0.005108635 0.0736452 1983 0.22380984 0.18804431 0.25957537 1984 0.04799352 0.011540157 0.08444688 1985 0.049729742 0.015663432 0.08379605 1986 0.09568698 0.064408004 0.12696595 1987 0.2430264 0.21218552 0.27386728 1988 0.2821517 0.24703526 0.31726816 1989 0.1792503 0.14449841 0.21400218 1990 0.36058238 0.32456872 0.39659604 1991 0.33889654 0.30403617 0.3737569 1992 0.12489683 0.09088211 0.15891157 1993 0.16570719 0.12821673 0.20319764 1994 0.23354979 0.19841295 0.26868662 1995 0.37686613 0.34365577 0.41007653 1996 0.27668938 0.24318002 0.31019875 1997 0.4223085 0.39009085 0.4545262 1998 0.5773417 0.54306877 0.61161464 1999 0.32448497 0.29283476 0.3561352 2000 0.33108476 0.29822785 0.36394167 2001 0.48928034 0.4580683 0.5204924 2002 0.5434665 0.51278186 0.57415116 2003 0.54417014 0.5112426 0.5770977 2004 0.46737072 0.43433833 0.5004031 2005 0.6068625 0.5757053 0.63801974 2006 0.5725527 0.541973 0.60313237 2007 0.5917013 0.56135315 0.62204945 2008 0.46564984 0.43265733 0.49864236 2009 0.5967816 0.5652556 0.62830764 2010 0.6803714 0.6490759 0.7116668 2011 0.53769773 0.5060012 0.5693943 2012 0.57760704 0.54485524 0.61035883 2013 0.6235753 0.5884838 0.6586669 2014 0.67287165 0.63890487 0.7068384 2015 0.8251144 0.7912871 0.85894173 2016 0.9329271 0.9017635 0.96409065 2017 0.84517425 0.81477475 0.87557375 2018 0.76265407 0.73105204 0.79425603 2019 0.8910726 0.85678726 0.9253579 2020 0.9229205 0.8882981 0.95754296 2021 0.76190555 0.7254577 0.7983534 2022 0.8013053 0.76525766 0.837353 2023 1.0858556 1.0421426 1.1295687]; yearHadCRUT5A = M(:,1); HadCRUT5A = M(:,2); HadCRUT5AL = M(:,3); HadCRUT5AU = M(:,4); %add 2024 HadCRUT5 %simply averaging across months as we see no significant seasonality %unlike the CO2 %but will weight the 2024 value by n/12 %so only equal weight in the regression when the annual sample is full M = [202401 1.151554 1.1032058 1.1999023 202402 1.2902366 1.2416096 1.3388637 202403 1.251485 1.1997101 1.30326 202404 1.2052884 1.1572673 1.2533094 202405 1.0745293 1.0205593 1.1284993 202406 1.1154083 1.0582764 1.1725402 202407 1.1379998 1.0766809 1.1993186 202408 1.2395844 1.1742324 1.3049364 202409 1.1410455 1.0814147 1.2006761]; yearHadCRUT5A = [yearHadCRUT5A;2024]; HadCRUT5A = [HadCRUT5A;mean(M(:,2));]; HadCRUT5AL = [HadCRUT5AL;mean(M(:,3))]; HadCRUT5AU = [HadCRUT5AU;mean(M(:,4))]; %temp measurement uncertainties to be used as regression weights %could do better in future by digging into those HadCRUT5 ensembles HadCRUT5Arang = HadCRUT5AU-HadCRUT5AL; HadCRUT5Asig = HadCRUT5Arang/4; %95th percentile hence assuming +/- 2 sig. %1850-1900 baseline tempanomalybaseline = mean(HadCRUT5A(1:50)); tempanomalybaselinesig = std(HadCRUT5A(1:50)) + mean(HadCRUT5Asig(1:50)); %------------------------------- %ANALYSIS %------------------------------- %construct interpolated single pCO2 record yearCO2 = (min(yearLawDome):max(yearMaunaLoa))'; CO2 = zeros(length(yearCO2),1)*NaN; %fill in Law Dome data for i = 1:length(yearCO2) j = find(yearLawDome == yearCO2(i)); if isempty(j)==1 CO2(i,1) = NaN; elseif isempty(j)==0 CO2(i) = LawDomeCO2(j); end end %fill in Mauna Loa data i = find(yearCO2>=min(yearMaunaLoa)); CO2(i) = MaunaLoaCO2; %how many points are interpolated? i = find(yearCO2>=1850 & yearCO2<=1958); j = find(isnan(CO2(i))==1); R = length(j)/length(i); %interpolate missing CO2 data with smoothed interpolation a la Enting i = find(isnan(CO2)==0); [CO2i] = interp1(yearCO2(i),CO2(i),yearCO2,'pchip');%plot(yearCO2,CO2,'ko',yearCO2,CO2i,'r');return [nvr,opts,parse] = irwsmopt([CO2i],1,'ml'); [CO2si] = irwsm(CO2i,1,nvr/1e3); %approximates Enting's 20 year spline smoothing %retain pre-1958 data i = find(isnan(CO2)==0); CO2si(i) = CO2(i); %retain post 1958 data i = find(yearCO2>1958); CO2si(i) = CO2(i); %CO2 measurement uncertainties CO2sig = zeros(size(CO2)); i = find(yearCO2<1959); CO2sig(i) = 2/2; %assume Etheridge quoted uncertainty of +/- 2ppmv - assume +/- 2sig? i = find(yearCO2>=1959); CO2sig(i) = 0.12/2; %see quoted Mauna uncertainty of +/- 0.12ppmv - assume +/- 2sig? %observed CO2 baselines i = find(yearLawDome<=1700); CO2baselinemu = mean(LawDomeCO2(i)); CO2baselinesig = std(LawDomeCO2(i)) + 1; %add Etheridge quoted measurement uncertainty of +/- 2ppmv - assume +/- 2sig disp('pre-1700 CO2 baseline') disp([CO2baselinemu CO2baselinesig*2]) i = find(yearCO2>=1850 & yearCO2<=1900); CO2baselinemu1850_1900 = mean(CO2si(i)); nvr = irwsmopt(CO2si(i),1); CO2baselinetrend = irwsm(CO2si(i),1,nvr); CO2baselinesig1850_1900 = std(CO2si(i)) + 1 - std(CO2baselinetrend); %add Etheridge quoted measurement uncertainty of +/- 2ppmv - assume +/- 2sig but subtract trend variance disp('1850-1900 CO2 baseline') disp([CO2baselinemu1850_1900 CO2baselinesig1850_1900*2]) %difference between the pre-1700 and 1850-1900 CO2 baselines N = 1e4; R1 = randn(N,1)*CO2baselinesig+CO2baselinemu; R2 = randn(N,1)*CO2baselinesig1850_1900+CO2baselinemu1850_1900; CO2baselinediffmu = CO2baselinemu1850_1900-CO2baselinemu; CO2baselinediffsig = std(R2-R1); %print baselines disp('difference between1850-1900 CO2 baseline') disp([CO2baselinediffmu CO2baselinediffsig*2]) disp('1850-1900 HadCRUT5 baseline') disp([tempanomalybaseline tempanomalybaselinesig*2]) %---------------------------- % REGRESSION - 1850 - 2023.8 %1. define x, y and weights i = find(yearCO2>=1850); x = CO2i(i)-CO2baselinemu; %we could feed the baseline uncertainty into the regression here y = HadCRUT5A; wy = 1./(HadCRUT5Asig).^2; wy(end) = wy(end)*10/12; %weight last year by completness of sample %2. iterative WLS A = [1 0 0]; for i = 1:4 xf = filter(A,1,x); yf = filter(A,1,y); [Xf,STDX,MSE,Sf] = lscov([xf ones(size(x))],yf,wy); m = Xf(1); c = Xf(2)/sum(A); Xfsig = sqrt(diag(Sf)); msig = Xfsig(1); csig = Xfsig(2)/sum(A); yhat = m*x + c; e = y-yhat; %th = ar(e,1);A = polydata(th); th = mar(e,1); [A,B,C,P] = getpar(th); end %finalised 1850-2023.8 estimates of senstivity and offset AA = A; PP = P; MM = m;MMsig = msig; CC = c;CCsig = csig; Sf = Sf; disp('1850-2023.8 regression result') disp([[m msig*2]*100;c csig*2]) AAsig = sqrt(diag(PP)); disp('AR(1) result') disp([AA(2)' AAsig*2]) %test residuals for normality ef = filter(AA,1,e); %decorrelate efsig = std(ef); [H,P] = adtest(ef,'alpha',0.001); disp('adtest result') disp([H P]) %test residuals for unit root [h,pValue,stat,cValue,reg] = adftest(e); %disp('e is not I(1)-which you know given ar(1)<<1') %---------------------------- % REGRESSION - 1850 - 2000 %1. define x, y and weights i = find(yearCO2>=1850 & yearCO2<=2005); x = CO2i(i)-CO2baselinemu; %we could feed the baseline uncertainty into the regression here i = find(yearHadCRUT5A>=1850 & yearHadCRUT5A<=2005); y = HadCRUT5A(i); wy = 1./(HadCRUT5Asig(i)).^2; %2. iterative WLS A = [1 0 0]; for i = 1:4 xf = filter(A,1,x); yf = filter(A,1,y); [Xf,STDX,MSE,Sfs] = lscov([xf ones(size(x))],yf,wy); m = Xf(1); c = Xf(2)/sum(A); Xfsig = sqrt(diag(Sfs)); msig = Xfsig(1); csig = Xfsig(2)/sum(A); yhat = m*x + c; e = y-yhat; %th = ar(e,1);A = polydata(th); th = mar(e,1); [A,B,C,P] = getpar(th); end %finalised 1850-2000 estimates of senstivity and offset AAs = A; PPs = P; MMs = m;MMssig = msig; CCs = c;CCssig = csig; Sfs = Sfs; %----------------------------- %simulate full stochastic model explicitly as don't know the analytical %solution for CIs using weights AND stochastic baseline %adding temp. weight effects on post MCS %start from 1700 j = find(yearCO2>=1850); yearMCS = yearCO2(j); M = length(j); N = 1e5; xMCS = zeros(M,N); eMCS = zeros(M,N); GMSTAMCS = zeros(M,N); HIW1MCS = zeros(M,N); HIW2MCS = zeros(M,N); SOMCS = zeros(M,N); %sample parameter covariance T1 = chol(Sf); N1 = randn(2,N); pse = (T1'*N1)'; aMCS = AA(2)+randn(N,1)*AAsig; mMCS = MM + pse(:,1); cMCS = CC + pse(:,2); %could include AR(1) uncertainty %MCS iterate for i = 1:N %check if CO2 baseline uncertainty active xMCS(:,i) = CO2i(j) - CO2baselinemu + randn*CO2baselinesig + randn(size(j)).*CO2sig(j); eMCS(:,i) = filter(1,[1 aMCS(i)],randn(M,1)*efsig).*HadCRUT5Asig./mean(HadCRUT5Asig); %AR(1) stochasticity weighted by normalised HadCRUT5 uncertainty GMSTAMCS(:,i) = mMCS(i)* xMCS(:,i) + cMCS(i) + eMCS(:,i); HIW1MCS(:,i) = mMCS(i)* xMCS(:,i) + cMCS(i) - CC;% + eMCS(:,i); %pre-1700 HIW trend HIW2MCS(:,i) = mMCS(i)* xMCS(:,i) + cMCS(i) - CC - (tempanomalybaseline + randn*tempanomalybaselinesig - CC) + eMCS(:,i); %1850-1900 HIW trend SOMCS(:,i) = (HadCRUT5A-CCs+randn*CCssig + randn(size(HadCRUT5Asig)).*HadCRUT5Asig)./xMCS(:,i); end %percentiles %for HIW and HICO2 use parametric mean +/-2 sig for consistency with regression results %otherwise use non-parametric percentiles through sorting %differences are otherwise small HICO250 = prctile(xMCS',50)';[HICO250s,i]=sort(HICO250); HIW105 = prctile(HIW1MCS',05)';HIW105s = HIW105(i); HIW150 = prctile(HIW1MCS',50)'; HIW195 = prctile(HIW1MCS',95)';HIW195s = HIW195(i); HIW205 = prctile(HIW2MCS',05)'; HIW250 = prctile(HIW2MCS',50)'; HIW295 = prctile(HIW2MCS',95)'; GMSTA05 = prctile(GMSTAMCS',05)';GMSTA05s = GMSTA05(i); GMSTA50 = prctile(GMSTAMCS',50)'; GMSTA95 = prctile(GMSTAMCS',95)';GMSTA95s = GMSTA95(i); SO05 = prctile(SOMCS',05)'; SO50 = prctile(SOMCS',50)'; SO95 = prctile(SOMCS',95)'; e05 = prctile(eMCS',05)'; e50 = prctile(eMCS',50)'; e95 = prctile(eMCS',95)'; %COL = [186 124 124]/255; COL = [175 100 100]/255; %-------------------------------------------------------------------------- %CO2-temperature plot xmin = 250;xmax = 450;xinc = 50; ymin = -.5;ymax = 2;yinc = 0.5; figure(1);clf hold on i = find(yearCO2<=1958 & yearCO2>=1850); j = find(yearHadCRUT5A<=1958); plot(CO2i(i),HadCRUT5A(j)-CC,'k.','markersize',10) i = find(yearCO2>1958); j = find(yearHadCRUT5A>1958); plot(CO2i(i),HadCRUT5A(j)-CC,'k+','markersize',4) h = patch([HICO250s;flipud(HICO250s)]+CO2baselinemu,[GMSTA05s;flipud(GMSTA95s)]-CC,[1 1 1]);set(h,'facecolor',[1 1 1]*.5,'edgecolor',[1 1 1]*.5,'facealpha',0.35,'edgealpha',0.35); %h = patch([HICO205;flipud(HICO295)]+CO2baselinemu,[GMSTA05;flipud(GMSTA95)]-CC,[1 1 1]);set(h,'facecolor',[1 1 1]*.5,'edgecolor',[1 1 1]*.5,'facealpha',0.35,'edgealpha',0.35); plot([CO2baselinemu;CO2i(end)],[0;MM*(CO2i(end)-CO2baselinemu)],'color',COL,'linewidth',1) plot([CO2baselinemu CO2baselinemu],[ymin 0],'k--') plot([xmin xmax],[1.5 1.5],'k--') plot([xmin CO2baselinemu],[0 0],'k--') plot([HICO250s(end) HICO250s(end)]+CO2baselinemu,[ymin 1.5],'k--') hold off flab = gca; flab.XMinorTick = 'on'; flab.YMinorTick = 'on'; xtickformat('%.0f'); ytickformat('%.1f'); %flab.XTick = [xmin:xinc:xmax]; %flab.XTickLabel = [xmin:xinc:xmax]; %flab.YTick = [1 10 100 1000 10000]; %flab.YTickLabel = '';%[10 10 10 10 10]; flab.FontSize = 20; flab.FontName = 'corbel light'; flab.FontWeight = 'normal'; xlabel('CO_2 (ppm)','fontsize',20,'fontname','corbel light') ylabel('warming ({\circ}C)','fontsize',20,'fontname','corbel light') box('on') axis([xmin xmax ymin ymax]) axis('square') %------------------------------ %Sensitivity change i = find(yearMCS>=2005); xmin = 2005-2;xmax = 2024+2;xinc = 2; ymin = 0.5;ymax = 1.5;yinc = 0.1; figure(2);clf hold on %plot(yearCO2(i),SO(i,2)*100,'col',COL*.8) h = patch([yearMCS(i);flipud(yearMCS(i))],[SO05(i);flipud(SO95(i))]*100,[1 1 1]); set(h,'facecolor',[1 0 0]*.5,'edgecolor',[1 0 0]*.5,'facealpha',0.25,'edgealpha',0.25); plot([xmin+2 xmax-2],[MMs-2*MMssig MMs-2*MMssig]*100,'--','col',COL*.8) plot([xmin+2 xmax-2],[MMs+2*MMssig MMs+2*MMssig]*100,'--','col',COL*.8) plot([2005 2005],[ymin ymax],'k--') plot([2024 2024],[ymin ymax],'k--') hold off flab = gca; flab.XMinorTick = 'off'; flab.YMinorTick = 'on'; xtickformat('%.0f'); ytickformat('%.2f'); flab.FontSize = 20; flab.FontName = 'corbel light'; flab.FontWeight = 'normal'; box('on') axis([xmin xmax ymin ymax]) axis('square') xlabel('year','fontsize',20,'fontname','corbel light') ylabel('sensitivity ({\circ}C per 100 ppm)','fontsize',20,'fontname','corbel light') %------------------------------ %2024 estimate %2024 HICO2 disp('2024 HICO2') disp([HICO250(end) HICO205(end) HICO295(end)]) %2024 HIW % i = find(yearCO2>=1700); % HICO2mu = CO2i(i) - CO2baselinemu; % HIW1mu = MM*HICO2mu; %i.e. the parametric estimate for consistency % HIW2mu = MM*HICO2mu - (tempanomalybaseline - CC); %i.e. the parametric estimate for consistency % disp('2024 R-BHIW') % disp([HIW1mu(end) HIW105(end) HIW195(end)]) % disp('2024 HIW using 1850-1900 as baseline') % disp([HIW2mu(end) HIW205(end) HIW295(end)]) %probability that 1.5 HIW is exceeded [HIW12024s,i]=sort(HIW1MCS(end,:)); p = (1:length(i))'./(1+length(i)); i = find(HIW12024s>=1.5); P15 = 1-p(i(1)); disp('P>1.5 degC for R-B pre-1700 HIW') disp(P15) [HIW22024s,i]=sort(HIW2MCS(end,:)); p = (1:length(i))'./(1+length(i)); i = find(HIW22024s>=1.5); P15 = 1-p(i(1)); disp('P>1.5 degC for R-B 1850-1900 HIW') disp(P15) [n,x] = hist(HIW1MCS(end,:),1e2); x = x'; p = 100*n'/sum(n); xmin = 1;xmax = 2;xinc = 0.1; ymin = 0;ymax = 5;yinc = 1; figure(3);clf hold on h = patch([x;flipud(x)],[p;zeros(size(p))],[1 1 1]); set(h,'facecolor',[1 0 0]*.5,'edgecolor',[1 0 0]*.5,'facealpha',0.5,'edgealpha',0.5); h = patch([x-(tempanomalybaseline-CC);flipud(x-(tempanomalybaseline-CC))],[p;zeros(size(p))],[1 1 1]); set(h,'facecolor',[1 0 0]*.5,'edgecolor',[1 0 0]*.5,'facealpha',0.25,'edgealpha',0.25); plot([1.5 1.5],[0 ymax],'k--') hold off flab = gca; flab.XMinorTick = 'off'; flab.YMinorTick = 'on'; xtickformat('%.1f'); ytickformat('%.0f'); flab.FontSize = 20; flab.FontName = 'corbel light'; flab.FontWeight = 'normal'; box('on') axis([xmin xmax ymin ymax]) axis('square') xlabel('warming ({\circ}C)','fontsize',20,'fontname','corbel light') ylabel('probability (%)','fontsize',20,'fontname','corbel light') return %--------------------END-----------------------