-1=e^ipi

But there is an overlap between instrumental data and the temperature reconstructions (ex. Mann et al. or Marcott et al.). Why can't one use the overlap to estimate the reconstructed data relative to the instrumental data? I can see chaos causing significant climatic variations on the order of decades (although most of that can be accounted for by certain factors such as AMO, PDO or ENSO). But being significant on the order of centuries? I don't buy it. If you want to offer me a physical mechanism as to why this would occur, I would be interested. Also, in my models I allow for error, so I just need sufficient data to overcome the uncertainty that is due to this chaos.

Update: My earlier claim that the AMO and the PDO can be modeled as 3 and 5 sinusoids respectively may not be correct. These oscillations may not have any periodic features at all; they may simply be random noise + autocorrelation and the observed cyclical behaviour is a result of this. If I take the AMO data, detrend a temperature response, a constant and autocorrelation, and then take the fourier transform of the residual, only the peak period of ~70 years remains. If I do the same thing with the PDO data, the PDO residual has no clear periodic tendencies. - - - - - I think I have a solution to a problem I was having earlier on how to detrend the El Nino Southern Oscillation (ENSO) from the temperature data set. The problem is that temperature change due to ENSO may be a lagged response. However the functional form of this response is unknown. If I assume a functional form, then I run the risk of unaccounted for specification error. Also if the functional form is difficult to estimate (even a simple exponential lag may be difficult) then it adds to the non-linearities of my model and makes it more difficult to have a convergent solution. Van Hateren assumed that temperature response to ENSO is exponentially lagged with a decay time of 4 months. Kevin Cowtan uses a simple lag of 4 months (http://www.ysbl.york.ac.uk/~cowtan/applets/nbox/nbox.html). This suggests that the time scale of response is much less than 1 year. So if I simply have both an ENSO term and an ENSO term lagged by 1 year in my model (or if I am using monthly data, maybe 12 ENSO terms with lags varying between 0 and 1 years) then this should give a good approximation of the temperature response to ENSO. Because of the a priori information and the fact that my degrees of freedom is reduced by adding in extra lagged terms, I shouldn't have significant unaccounted for specification error. It also enters the functional form linearly, which makes it easier to estimate. I can probably do the same thing for the temperature response to AMO or PDO since the natures of these responses are similar. Also, when I try to detrend temperature from AMO, PDO or MEI in order to avoid reverse causality, perhaps I should have lagged temperature terms (up to 1 year) to account for the possibility of a lagged temperature response. - - - - - I managed to construct a monthly data set from 1866-2012 and accounted for all of the issues associated with unequal days in each month + more accurate data at later times (plus seasonally detrending and removing temperature response from AMO, PDO and MEI). It didn't really improve the estimates by much (even after I varied the model a few times), probably because even though I do have more data points, the variation in the parameters in data is still small (solar irradiance, greenhouse gas forcing and SO2 are all strongly correlated over this period). The biggest uncertainty seems to be with how strong solar activity is relative to greenhouse gas forcing. Another issue is that if I use monthly data, convergence to a unique estimate is a lot harder (since my model is nonlinear) and it takes far more computational time. - - - - - I found reconstructions of 5 more long-lived greenhouse gases (CFC-11, CFC-12, CFC-113, CCl4, SF6) http://cdiac.ornl.gov/oceans/new_atmCFC.html. I can calculate their radiative forcing effect using https://www.ipcc.ch/publications_and_data/ar4/wg1/en/ch2s2-10-2.html. Including these factors should give a better estimate of greenhouse gas forcing than just using CO2, CH4 and NO2. In particular, since CFC concentrations have been decreasing since the 90's, this can help explain some of the slowdown in global warming over the past 2 decades.

The literature is in disagreement with itself. That's one of the reason's the IPCC's uncertainty range for climate sensitivity (1.5 C - 4.5 C) is so large. Instrumental climate sensitivity estimates (and to a lesser extent the paleo estimates) consistently give low climate sensitivity estimates (1.5 C - 3.0 C), where as the general circulation computer models are consistently giving high climate sensitivity estimates (3.0 C - 4.5 C). Of course, a lot of this perceived uncertainty is partly due to individuals not analyzing the data properly and making gross oversimplifications (the Loehle paper in the original post is a good example, but I would argue that even individuals like James Hansen do this). In addition, many people are underestimating the uncertainty of their estimates because their models may have a lot of specification error (this is especially true for GCMs). Also, the definition of 'equilibrium climate sensitivity' is somewhat ambiguous, which adds to the uncertainty of the estimates. I'm trying to do time series analysis of the empirical data properly to see if I can get a decent estimate of climate sensitivity (and the impulse response function) without making oversimplifications that give large specification error. Are the empirical estimates in the literature low because of specification error? Or are the GCM estimates in the literature high because of specification error? I don't know.

e^ix gives you a circle on the complex plane if you vary x from 0 to 2*pi.

@ TimG How about Marcott et al. is that reasonable to use or not? I can take into account the issues of the reconstructed temperatures having a lower time resolution and greater uncertainty (by simply doing a weighted regression that takes measurement error into account). But if the temperature estimates are biased or the uncertainty is underestimated then that is a concern. Two reasons why using reconstructed data may be useful are because it can give an idea of the long term temperature responses (that may occur due to a change in CO2 levels) and it can help me better distinguish between temperature changes due to solar activity and temperature changes due to greenhouse gases. If there is only 130 years of empirical data then it becomes difficult to estimate temperature changes over a longer time scale than this. Also, since solar activity, greenhouse gas forcing and SO2 levels are very correlated over the last 130 years, there isn't enough data (or data of good quality) to be able to give a good estimate of climate sensitivity if I allow the temperature response to solar activity or SO2 levels to vary freely with respect to the temperature response to changes in CO2 levels. I could impose restrictions such as Van Hateren did, but then I will have greater specification error and the estimate may be biased. Anyway, I don't mind giving the benefit of the doubt to the climate alarmists for the sake of argument. As far as I can tell, I don't see evidence of high climate sensitivity estimates from the empirical data. The instrumental data seems to support an equilibrium climate sensitivity between 1.5 C and 3 C. So I'm seeking alternative 'explanations' as to why my estimates are consistently low. Also, if Mann is biased downward, then that will mean that the estimate of the temperature response to a change in forcing will be low for longer response times, which may actually lead to an underestimate of climate sensitivity. So the climate alarmists can't have it both ways. If the temperature variation before industrial times is small, then this suggests that the long-term temperature response to a change in forcing is small (but at the same time would suggest that a greater percentage of the warming observed over the past 130 years is due to greenhouse gases).

Perhaps I should use monthly data from 1876-2012 instead of annual data. That might help reduce the uncertainty. In addition, perhaps I should weight months by the inverse square of their standard errors (which I can get from http://www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/HadCRUT.4.3.0.0.monthly_ns_avg.txt). Another option is to use reconstructed temperature data from earlier points in time. Mann et al. 2008 seems to be the best reconstruction for 500-1850 AD, and Marcott et al. 2013 seems to be the best reconstruction over the Holocene (both are shown in the image below). If I want to get really crazy I can use Pleistocene reconstructed temperatures by Shakun and Carlson 2010. Note that the blue area is the 1 sigma uncertainty for Marcott et al., so don't put too much confidence in it. The Marcott et al. data has far less time resolution (20 years at best) and more uncertainty. Also, if I were to use data for the entire holocene then I would have to take into account Milankovitch cycles (which explain why Holocene temperatures have been decreasing for the past 8000 years). There are also more reconstructed data sets available for the Mann et al. time period than the Marcott et al. time period. Also, I want to retract an earlier claim I made about current global temperatures relative to the medieval warm period. Current global temperatures are warmer than the majority of the medieval warm period, but based on the uncertainty of Mann et al. chances are that at least 1 year in the medieval warm period was warmer than 2014.

So apparently human emitted SO2 has a ridiculously fast decay time of like half a day. http://www.sciencedirect.com/science/article/pii/0004698175901286 Therefore, human emitted SO2 forcing should be roughly proportional to SO2 emissions. Finding global SO2 emission data is fairly difficult. However, I managed to find historic SO2 emission data (http://sedac.ciesin.columbia.edu/data/set/haso2-anthro-sulfur-dioxide-emissions-1850-2005-v2-86/data-download) and recent SO2 emission data (http://iopscience.iop.org/1748-9326/8/1/014003/media/erl441620suppdata.pdf). Both data sets are constructed by the same people (Smith et al) although there is a slight discrepancy from 2000-2005 during the overlap. I can assume that the 2000-2011 data set is more accurate and simply increase the 2000-2011 data by 541 GgS/y to make them comparable. For a rough value of 2012 SO2 emissions, I can estimate this using the linear trend between 2010-2011. SO2 emissions are the only thing that can explain my low estimates of climate sensitivity (other than climate sensitivity actually being low) since they have a significant cooling effect. However, if I put SO2 as a climate variability index or if I use SO2 directly as a forcing (and either treat it as a free parameter or fix it to something reasonable such as -0.005 W/m^2 per TgS/y http://onlinelibrary.wiley.com/doi/10.1029/2007JD008683/pdf) I get a higher climate sensitivity, but the uncertainty becomes ridiculously large (example: 95% confidence interval of [0.05,184.93]C with best estimate of 4.25C). It looks like I'll either need more data or better model specification if I want a reasonable estimate of climate sensitivity. - - - - - One more thing, I was looking over the Van Hateren 2012 paper again and I think it both underestimates climate sensitivity and underestimates the uncertainty for a number of reasons: - To get the relative effect of changes in solar irradiance in W/m^2 and greenhouse gas forcing in W/m^2, Van Hateren takes 1 - 0.3 (which he/she assumes to be the albedo of the earth) and divides by 4 (which is the ratio of the surface area of a sphere to the area of a circle of equivalent radius). However, as I mentioned earlier, changes in solar irradiance affect equatorial regions more due to receiving more direct sunlight. And due to the Stefan-Boltzman law, the temperature change due to changes in solar irradiance should be less than what Van Hateren is assuming (I explained this in more detail in post #269, where I gave a lower bound of 0.0441). Not to mention 0.3 for the Earth's albedo is likely an underestimate. Although there is also the effect of cosmic rays. Overall, this assumption about the relative strength of solar activity vs greenhouse gas activity should cause an underestimate of climate sensitivity. - Van Hateren allows for the coefficients of his impulse response function to be negative, which as I explained earlier can give a nonsense impulse response function that underestimates climate sensitivity. Although this doesn't seem to be an issue by the results of figure 6 in Van Hateren's paper. - Van Hateren detrends data with the Multivariate ENSO index (MEI) but doesn't take into account the issue of reverse causality. If global warming causes MEI to increase, then this means that Van Hateren is underestimating climate sensitivity. - There is a lot of model specification error in what Van Hateren does, which is not taken into account. I'm not going to go into all the details, but some of the stuff I have done earlier in this thread is comparable with what Van Hateren has done and if I relax my assumptions to try to take into account specification error, my uncertainty increases by quite a bit. Also, there are a lot of coefficients that Van Hateren assumes to be true without much justification that he uses in detrending and in other things (and the uncertainty of these coefficients is not taken into account).

Looking at the variability indices, I think that I will go with this for a while: - Use a 3 sinusoidal model to detrend AMO from temperature. - Use a 5 sinusoidal model to detrend PDO from temperature. - Use SOI instead of MEI, since ENSO is too chaotic to predict and SOI should be temperature detrended. Perhaps I need to determine the optimal number of sinusoids using Akaike's Information Criterion. Doing this gives me a 95% confidence interval of [0.71, 3.79]C with a best estimate of 2.00 C. This is a slight reduction in uncertainty from the earlier estimates. However, the residual from this estimate still shows some features of the above 3 climate oscillations. One possibility is that there is a lagged response in global temperatures to these climate oscillations. One solution is to impose an exponential lagged effect of the climate oscillations on global temperatures. Using a decay time of half a year gives a 95% confidence interval of [0.80, 2.97]C with a best estimate of 1.79 C. This gives an even smaller confidence interval. Perhaps I should have this decay time as a free parameter in the model, or I should treat changes in these climate oscillations as forcings. But having a lagged effect of climate indices seems to be a good idea. One possible explanation for my low estimates is that I'm not fully taking human emitted aerosols into account. So I will need to look into this.

I thought I would try to fit simple models to AMO, PDO and MEI (to get better detrended data, and maybe to use as projections). So I took AMO, PDO and MEI and looked at the fourier transform. AMO has strong peaks around ~70 years, 24 years and 10 years. PDO has strong peaks around 40 years, 27 years, 13 years and 5.7 years. MEI has ~70 years, 10 years, 7 years, 4 years and 3 years (as well as others). I tried using the gauss newton algorithm to fit these indices to a constant + a linear factor + temperature + a few sinusoids. The AMO can be approximated very well by 2-3 sinusoids. The PDO can be approximated okay by 3 (but it might need more; My computer only has so much computational power, so I haven't done 4 yet). MEI is very chaotic so is difficult to predict without a very large number of sinusoids.

I thought I would try changing my variability indices to see if it would significantly change the estimate. The problem with using Multivariate ENSO Index (MEI), Pacific Decadal Oscillation Index (PDO) and Atlantic Multidecadal Oscillation Index (AMO) is that they all use temperatures, so there is the issue that global warming may change these indices (reverse causality issue). To try to prevent this, I have take each of these indices are regress it on a constant, on a linear term and on temperature (then I take the residual as my detrended index). Using the detrended PDO changes my 95% confidence interval to [0.65,4.67]C with a best estimate of 1.74 C. Replacing the SOI with the detrended MEI changes my 95% confidence interval to [0.70,3.91]C with best estimate of 2.05 C. Replacing the NAO with the detrended AMO changes my 95% confidence interval to [0.66,4.26]C with best estimate of 2.15 C. My confidence interval hasn't changed much by modifying these indices. Maybe I need to construct models for how each of these indices vary over time to better detrend them (this could also be used to make temperature projections). One thing I should mention. For the original post, I didn't take into account the fact that there are other greenhouse gases than just CO2. about 76% of the change in greenhouse gas radiative forcing between 1950 and now is due to CO2. This suggests that my estimate of climate sensitivity from the original post should be 2.25 C rather than 2.95 C. Of course this still has the assumption of constant natural warming (which means that the estimate should be higher than 2.25 C).

Update: I've been reading a bit on the literature of the climate impulse response function. Seems like Fortunat Joos is one of the most cited people on this. This paper http://www.climate.unibe.ch/~joos/papers/hooss01cd.pdf is fairly interesting because it determines impulse response functions for many things (temperature change due to CO2, sea level rise due to CO2, CO2 decay rate towards equilibrium, etc.). On page 7 it gives a two exponential impulse response function for temperature response to CO2. One has a decay time of 12 years, and another with a decay time of 400 years. If one interprets that first decay time to correspond to the equilibrium climate sensitivity and the second decay time to correspond to the earth system sensitivity, then this suggests that earth system sensitivity is approximately 1/0.71 = 1.408 times equilibrium climate sensitivity. This is completely consistent with the claim by Lunt et al. 2012 (which used Pliocene paleoclimate data) that the ESS is approximately 1.4 times the ECS. If the above paper is correct about a 12 year decay corresponding well to equilibrium climate sensitivity, then if I use a 12 year decay, I should be able to get a good estimate of equilibrium climate sensitivity. Doing this gives me (1.74 +/- 0.18) C, this is pretty low (of course this uncertainty doesn't take into account the large specification error of choosing a 12 year decay time). This webpage (http://unfccc.int/resource/brazil/climate.html) gave a similar 2 decay time model, but with 8.4 year and 410 year decay times. Another thing I did was I took the impulse response function from the end of this paper (http://www.climate.unibe.ch/~joos/IRF_Intercomparison/Protocol_CO2_impulse_response_modelcomparison_v1.0.pdf) using a freeware program called GetData. I wanted to see what functional forms can represent this impulse response function. The simple exponential decay (1 - exp(-At)) is terrible as expected. However, (1 - exp(-A*t^B ) and (1 - exp(-(A - B*exp(-Ct))*t)) are decent. Maybe I could try using these rather than the Van Hateren approach, although the additional non-linearity will be a pain to deal with. Alternatively, perhaps I should just take one of these GCM impulse response functions as is, and see what climate sensitivity I get from the data under the assumption of this impulse response function. - - - - - So most of my recent estimates have been suggesting a climate sensitivity lower than 3 C. To play a climate alarmist advocate, one possible reason for this is that my regression model results in solar and volcano coefficients that may be too high (therefore, it overstates their importance relative to CO2 and gives a low estimate of climate sensitivity). Therefore, it might be interesting to choose lower bounds on what these two coefficients should be, and then try to estimate climate sensitivity from the data using a regression. For volcanic aerosols, as I explained earlier, the NASA website suggests that one can convert the volcanic aerosol data to forcing in W/m^2 by multiplying by a factor of -23. However, Van Hateren mentioned in his paper that the literature suggests that a forcing of volcanic aerosols only has about 50-60% the effect of an equivalent forcing of greenhouse gases. To play the climate alarmist advocate, I'll choose the lower of these values, which means my volcano coefficient should be no less than -11.5. Trying to put a lower bound on the solar constant is more tricky. Solar irradiance is correlated with cosmic rays, which may amplify the measured effect of solar irradiance, but since the sun's solar irradiance has a stronger effect in equatorial regions than polar regions, this unequal heating affect should result in a smaller temperature response than if the change in solar irradiance was distributed uniformly (due to the Stefan-Boltzman Law). As I argued earlier, the fact that the ratio of the surface area of a sphere to the area of a circle of equal radius is 4 and the fact that the earth has an albedo of 0.31 (earlier I said 0.35, but 0.31 seems to be a better value http://www.climatedata.info/Forcing/Forcing/albedo.html) suggests that the solar coefficient should be about 0.1725. To play the climate alarmist advocate, I'll pretend cosmic rays don't exist and try to calculate the affect of unequal heating to get a lower bound for the solar coefficient. The average direct solar irradiance directed upon a square meter of Earth at latitude φ is proportional to: ϑ*Integral(μ = 0 to 2π; [cos (φ) cos (sin-1 (cos (μ) sin (ψ0 ) ) ) sin (ψ) - sin (φ) cos (μ) sin ( ψ0) ψ ]dμ) Where ϑ = 1366 W/m^2 is the solar irradiance at 1AU from the sun, ψ0 = 23.4*π/2 is the Earth’s tilt in radians, and ψ = cos-1 (min(0,max(1,tan (φ) tan (sin-1 (cos (μ) sin ( ψ0) ))) ). I’ll spare you guys the derivation (since you guys complain my posts are too long). I’m missing a normalization constant here. In any case, I can normalize this and then make sure the average direct solar irradiance equals ϑ/4/π. Now consider a simplified representation of the earth, where the temperature distribution depends only on the latitude temperature profile and the rate of heat transfer between two adjacent square meters of the surface is proportional to the temperature difference. Then in equilibrium, one has: (1-α)*S(φ) + B = G*σ*T4(φ) + k*(d2T(φ)/dφ2 - tan(φ)*dT(φ)/dφ) Where α = 0.31 is the albedo of earth, S(φ) is the solar irradiance at latitude φ, B = 0.087W/m^2 is the heat flux due to the Earth’s internal energy, G is a factor due to greenhouse gases, σ is the Stefan-Boltzmann constant, k is a constant that determines the rate of heat transfer, and T is the temperature. Again, I’ll spare you guys the derivation. Now the above equations have no analytical solution. Though it can be approximated numerically with a computer. I need to determine k and G before I can see how the unequal heating of increasing solar irradiance compares to the equal heating of increasing greenhouse gases in terms of equilibrium temperature response. For an earth with approximately uniform temperature, integrating the second equation gives: (1-α)*ϑ/4/π + B = G*σ*T4 For a temperature of 288 K, this suggests that G = 0.1925. If I look at the temperature profile of Earth: http://www.roperld.com/science/graphics/TempVSLatitude.jpg At a latitude of 40 degrees, the temperature is approximately 286K, the first derivative of temperature with respect to latitude is approximately -0.7 C per degree, and the second derivative is approximately -0.1 C per degree squared. Putting this into the second equation gives: (1-α)*S(40*π/180) + B = G*σ*(286 K)4 + k*((-0.1K + 0.7K*tan(40*π/180))*(180/π)^2) This gives k = -8.78 * 10^-4 W/m2/K. The other thing I need is an initial temperature profile. I’ll start with a simple sinusoidal temperature profile of 273 + 26*cos(2φ) K. From here, I can write a computer program to start with this initial temperature profile and slowly change the temperature profile based on the imbalance of equation 2 until equilibrium is reached. I can then allow for k and G to vary and impose restrictions that the average temperature is 288K and the equatorial temperature is 300 K to get estimates of k and G. I get that k = -0.0646 and G = 0.1956. Now to see the equilibrium temperature response to a change in solar irradiance, I simply increase solar irradiance slightly (say by 0.1 W/m^2) and calculate the change in average temperature once the new equilibrium is reached. To see the equilibrium temperature response to a change in greenhouse gases, I simply decrease G such that the average of G*σ*T4 decreases slightly (say by 0.1 W/m^2) and calculate the change in average temperature once the new equilibrium is reached. I find that the effect of increasing solar irradiance is about 4.41% as effective as increasing greenhouse gas forcing by an equivalent W/m^2. So Solar = 0.0441 and Aerosol = -11.5 are my lower bounds for these coefficients. If I use this then the regression becomes: T(Y) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dFi*ρs(i,i) + ... + dF1700*ρs(i,1700)) + β0 + β1*LODY + β2*SOIY + β3*PDOY + β4NAOY + errorY Where dF is the sum of the forcing changes (dGHG + 0.0441*dS – 11.5*dA). This is fortunately a linear model, so is very easy to estimate. This model gives me a climate sensitivity of 0.97 +/- 0.59. The reason why I get this nonsense result is that some of the exponentials are negative. So I have to impose restrictions that they are positive. Fortunately, I have found a methodology that does this (thanks to stata's website) and allows me to appropriately measure the error. All I have to do is replace each γ with exp(ln(γ)) and then use the ln(γ) as my unknown coefficients. The regression equation becomes: T(Y) = Sum(s = 1 to 5; exp(lnγs)*((Sum i = 1876 to Y-1; dFi*ρs(i,i) + ... + dF1700*ρs(i,1700)) + β0 + β1*LODY + β2*SOIY + β3*PDOY + β4NAOY + errorY Of course, this is now a non-linear regression. But I can estimate this using Gauss-Newton. One issue is that if one of my γ's is supposed to be zero, then ln(γ) will approach minus infinity, which will eventually cause matlab to divide by zero (due to rounding error) and give me nonsense results. Therefore, I have to modify my code slightly to tell matlab not to make any of the lnγ's smaller than say -12 (which will mean that the γ is basically zero). Doing this should not result in an overestimation of my uncertainty (which is good). Anyway, I chose decay times of 0.5, 1.5, 4.5, 13.5, and 40.5 years (because I was having more computer-rounding issues with longer decay times). Note that my estimates of ln(γ) will all be normal, which means that each of my γ will be log-normal. This means that my estimate of climate sensitivity will be a sum of log-normals that are all correlated with each other (which I'll assume is roughly log-normal as well). This gives me a 95% confidence interval of [0.69,4.26] C for equilibrium climate sensitivity and a best estimate of 1.71 C. Unlike earlier results, I think that the confidence interval is well approximated here and there isn't much specification error on this model. Interestingly, this excludes the absurdly low climate sensitivity estimates by people like Christopher Monckton and excludes the very high end of the IPCC's confidence interval (which is 4.5C). Also, my residual is very correlated with MEI and AMO, so I might be able to reduce my uncertainty by a fair amount. One more thing, I forgot to multiply most of my confidence intervals by a factor of 5.35*ln(2) in the last 1-2 pages. Sorry about this mistake.

There is a potential issue in the model of post #260. If there is a significant amount of error for T(1876) - LOD(1876) - SOI(1876) - PDO(1876) - NAO(1876), then this will affect all of the results by a constant factor (which the regression will try to correct for with a weird impulse response function). I can modify the model by introducing a constant. The modified model is given below. Note that even though I increase the number of coefficients that need to be estimated by 1, the number of data points available to me also increases by 1, so the uncertainty of my estimates should not increase due to this change. T(Y) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1700*ρs(i,1700)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1700*ρs(i,1700)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1850*ρs(i,1850)))) + β0 + β1*LODY + β2*SOIY + β3*PDOY + β4NAOY + errorY I also found an error in my code where I was using solar irradiance data before 1876 instead of volcanic aerosol data before 1876. I've also modified my code such that I can calculate climate sensitivity given any set of decay times. My climate sensitivity estimates are now more robust and more accurate, but the outcome is fairly sensitive to the choice of the decay times. Here are some example results: Choosing the decay times to be 0.5, 2, 8, 32 years gives a climate sensitivity of (2.24 +/-0.34) C. Choosing the decay times to be 0.5, 1.5, 4.5, 13.5 and 40.5 years (factors of 3 difference) gives a climate sensitivity of (3.27 +/- 0.97) C. And since you guys kept saying 42, a decay time of 42 years gives a climate sensitivity of (3.28 +/- 0.17) C. It seems that a decay time of 128 years is too long to detect however. I would say that these confidence intervals are small enough such that non-linearities do not significantly affect them. Although, these confidence intervals do not account for the specification error in the choice of decay times. If anyone has a set of decay times they want me to test, just ask me to in this thread. Perhaps I should use some more a priori justified decay times as Bonam has suggested, or change the model to allow for the decay times to vary. Edit: actually these confidence intervals are off by a factor of 5.35*ln(2).

Nevermind, the 2.50 C estimate isn't valid. I checked the residual and it has a clear downward trend. My restricted estimation method could not deal effectively with one of the gammas being zero. In any case, the restriction suggested that the gammas for the decay times of 8 years and 128 years should be zero (or negative). I retried the model with only 3 decay times. This gave a climate sensitivity estimate of (1.28 +/- 0.03) C (again, I don't think the premises used to obtain the confidence interval are valid). Checking the residual, it seems to be correlated with the Multivariate ENSO Index (0.21 correlation coefficient) and the Atlantic Multidecadal Oscillation Index (0.42 correlation coefficient). So perhaps my current internal variability indices are not sufficient. The problem with using MEI or AMO is that these indices use temperature data and might be correlated with global warming. There are a few people that have tried to detrend this effect though. Sorry that I haven't been able to get something conclusive yet. Although, the closest results to conclusive that I have are suggesting a lower climate sensitivity.

I think I may have confused the North Atlantic Oscillation with the Atlantic Multidecadal Oscillation...

So I tried increasing the number of decay times. I now have 9: 0.5, 1, 2, 4, 8, 16, 32, 64, and 128 year decay times. This lead to a sensitivity of 0 +/- 9 C (there might be too many parameters to estimate). I tried to impose the restriction that the γ cannot be negative. This lead to a climate sensitivity of 3 C, but with a very large uncertainty (and my current method to evaluate uncertainty is probably invalid now that I've added the restrictions). I tried the restriction in the simpler model with only 5 decay times. This gives a climate sensitivity of (2.50 +/- 1.67) C. But again, the method to evaluate uncertainty is probably invalid. Edit: actually, the 2.50 C is invalid.

My impulse response function has a weird behavior since 2 of my γ are negative. The impulse response function is peaking after ~3 years. I should either make the set of γ more dense (such as have them spaced out by a factor of 2 rather than 4) or restrict my γ such that they cannot be negative when performing my estimation.

As a robustness check, I decided to remove the 4 internal variability variables and test: T(Y) - T(1876) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1700*ρs(i,1700)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1700*ρs(i,1700)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1850*ρs(i,1850)))) + errorY I get a climate sensitivity of (1.14 +/- 0.33) C.

Okay, I tested the model T(Y) - T(1876) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1700*ρs(i,1700)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1700*ρs(i,1700)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1850*ρs(i,1850)))) + β1*(LODY – LOD1876) + β2*(SOIY – SOI1876) + β3*(PDOY – PDO1876) + β4*(NAOY – NAO1876)) + errorY I get a climate sensitivity of (1.11 +/- 0.20) C. Seems low... But my estimates of the other coefficients seem reasonable.

Sorry, I've been writing NO2 instead of N2O the past two pages. Also, I found a discontinuity in my methane data, which I have corrected for. I've thought of a way to resolve the issue I've been having with the dG. The problem is that in 1876, it would not be reasonable to assume that global temperatures are in equilibrium since greenhouse gases and solar irradiance were steadily increasing since the 1700s. However, it might be reasonable to assume that global temperatures were in equilibrium in 1700 (at the end of the Maunder Minimum) and that any changes in radiative forcing that occurred before 1700 has negligible impact on changes in global temperatures after 1876. If this is the case then I can rewrite my regression equation as: T(Y) - T(1876) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1700*ρs(i,1700)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1700*ρs(i,1700)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1700*ρs(i,1700)))) + β1*(LODY – LOD1876) + β2*(SOIY – SOI1876) + β3*(PDOY – PDO1876) + β4*(NAOY – NAO1876)) + errorY Which means that the only thing left to do to estimate the above model is to obtain radiative forcing data from 1700-1875. Unfortunately, the volcanic aerosol data only goes back to 1850. However, this shouldn't be a big deal because volcanic aerosols quickly decay and volcanic activity there are no long term trends in volcanic activity, so volcanic activity prior to 1850 (even the 1815 Mount Tambora eruption that lead to the year without a summer) should have negligible impact on temperature changes after 1876. So I'll exclude changes in volcanic aerosols before 1850, and my regression equation becomes: T(Y) - T(1876) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1700*ρs(i,1700)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1700*ρs(i,1700)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1850*ρs(i,1850)))) + β1*(LODY – LOD1876) + β2*(SOIY – SOI1876) + β3*(PDOY – PDO1876) + β4*(NAOY – NAO1876)) + errorY The solar irradiance data I was using (http://lasp.colorado.edu/data/sorce/tsi_data/TSI_TIM_Reconstruction.txt) extends back to 1610, so I can just use the same data. The law dome CO2 data set that I was using (http://cdiac.ornl.gov/ftp/trends/co2/lawdome.combined.dat) goes back to 948. One problem is that the 20 year smoothed spline data that I was using only goes back to 1832. I can use the 75 year smoothed spline data before then, but then there is a discontinuity in my data set. From 1850-1870, the 75 smoothed spline data is on average 0.362857143 ppm more than the 20 year smoothed spline data. Therefore, I increase the 75 year smoothed spline data by 0.362857143 ppm to make it comparable. I then use the modified 75 year smoothed spline from 1700-1850 and the 20 year smoothed spline from 1851-1958. I then increase the ice-core data by 0.2765 ppm to make them comparable to the Mauna Loa data. For methane, the 75 year smoothed spline ice core data that I was using (http://cdiac.ornl.gov/ftp/trends/atm_meth/EthCH498B.txt) goes back until 1010, so I can extend my data back to 1700 using the same methodology I used earlier. For N2O, it is a bit tricky. The snow-pack data (ftp://daac.ornl.gov/data/global_climate/global_N_cycle/data/global_N_perturbations.txt) that I was using starts in 1756 and it is quite noisy prior to 1876. Ice core data for N20 is is available (ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/law/law2006.txt), but is still quite noisy and seems to give values lower than the snow-pack data. I take the ice-core data points from 1681 to 1900 and perform a quartic fit. If I compare the fitted values of the ice-pack data with the fitted values of the snow-pack data I find that in 1876, the snow-pack data is approximately 3.787136171 ppb higher than the ice-pack data. Therefore, I increase all of the icepack data by 3.787136171 ppb, use ice-pack data from 1700-1875, and use snow-pack data from 1876-1999. I then subtract 0.01283 ppb from both the snow-pack and ice-pack data to make it comparable to the 2000-2012 instrumental data. From this I can get estimates of greenhouse gas forcing from these 3 gases from 1700-2012.

So for the annual data from 1876-2012, the regression equation is: T(Y) - T(1876) = Sum(s = 1 to 5; γs*((Sum i = 1876 to Y-1; dGHGi*ρs(i,i) + ... + dGHG1876*ρs(i,1876)) + Solar*(Sum i = 1876 to Y-1; dSi*ρs(i,i) + ... + dS1876*ρs(i,1876)) + Volcano*(Sum i = 1876 to Y-1; dVi*ρs(i,i) + ... + dV1876*ρs(i,1876)) + dGs*(Sum i = 1876 to Y-1; ρs(i,1875)))) + β1*(LODY – LOD1876) + β2*(SOIY – SOI1876) + β3*(PDOY – PDO1876) + β4*(NAOY – NAO1876)) where Y is the year. ρ(i,j) = 0 if j > i ρ(i,j) = 1 + τ*(exp(-1/τ) –1) if i = j ρ(i,j) = τ*exp(-i/τ)*(exp(-1/τ) + exp(1/τ) - 2) if j < i Where I tried to estimate this using the technique of a sequence of linear regressions with alternating restrictions, it does not converge. When I try to perform the normal Gauss-Newton technique, it does not converge. When I modify the Gauss-Newton technique such that the change in the estimate between consecutive iterations is only 1% of the original Gauss-Newton technique, I do get initial convergence. I get that equilibrium climate sensitivity is 2.64 C. Note that in my model, equilibrium climate sensitivity is simply (γ1 + γ2 + γ3 + γ4 + γ5)*5.35*ln(2). However, if I tell matlab to keep going, I eventually get divergence and head to a case where sensitivity approaches zero and the dG terms get large. So perhaps I need to find a reasonable way to place restrictions on dG. As for the estimate of the error for climate sensitivity, it is ridiculously large, but I think that is mostly because my method of estimating error completely fails since the model is so non-linear that even the Gauss-Newton method does not converge properly. I think I need to read up more on convergence of non-linear regressions.

I found a mistake in my code which might have led to some of the nonsense results I was getting. I was using dμ = μ(endi, endi+1,dF) - μ(starti, endi,dF) instead of dμ = μ(endi, endi+1,dF) - T(0) when endi occurs at t=0 (see post #236 for more details). Edit: retried the regression of post #254. The results are basically the same.

Okay, I constructed a data set for 1876-2012. I removed atmospheric angular momentum from the model because there isn’t good AAM data prior to the 50s. However, AAM is strongly correlated with LOD and SOI, so it should not be a big deal. Temperature: http://www.cru.uea.ac.uk/cru/data/temperature/HadCRUT4-gl.dat CO2: ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_annmean_mlo.txt CH4: I use Cape Grim, Tasmania, Australia Instrumental data from 1985-2013 ftp://ftp.cmdl.noaa.gov/data/trace_gases/ch4/flask/surface/ch4_cgo_surface-flask_1_ccgg_month.txt and 1870-1984 spline smoothed ice core data from Antarctica and Greenland http://cdiac.ornl.gov/ftp/trends/atm_meth/EthCH498B.txt. I interpolate the ice core data for missing years and I decrease the ice core data by 24.12 ppb so that it is continuous with respect to the Cape Grim data. N2O: For N2O, I use global methane data from 2000-2014 ftp://ftp.cmdl.noaa.gov/hats/n2o/insituGCs/CATS/global/insitu_global_N2O.txt and I use snowpack Antarctica data from 1870-2004. ftp://daac.ornl.gov/data/global_climate/global_N_cycle/data/global_N_perturbations.txt . The data prior to the 90’s is relatively inaccurate and has missing data points, so I perform a cubic fit to the data prior to 1996. The cubic fit is relatively good (R^2 is 0.98) and intersects the data in the 90s around 1993. Therefore, I use the cubic fit for the data prior to 1994. For the 5 overlapping years (2000-2004) between the two data sets, the snowpack data is on average 0.01283 ppb higher than the later data set. Therefore, I reduce the snowpack data by 0.01283 ppb to make the data sets comparable. TSI: http://lasp.colorado.edu/data/sorce/tsi_data/TSI_TIM_Reconstruction.txt Volcanic Aerosols: http://data.giss.nasa.gov/modelforce/strataer/tau.line_2012.12.txt I thought I would try to get estimates of volcanic aerosols for October, November and December for 2012 (so that I have 1 more data point to work with if I want to use annual data). I’ll use the annual volcanic aerosol data from 1850-2011 and I’ll use the volcanic eruptions of VEI > 3 since 1850. I’ll define a variable V which is equal to the number of volcanic eruptions of VEI 4 plus 10 times the number of volcanic eruptions of VEI 5 plus 100 times the number of volcanic eruptions of VEI 6 for that year. Then I’ll regress dAOD(t)/dt = A + B*V(t) – k*AOD(t), where the first term corresponds to aerosol emissions due to volcanic eruptions of VEI less than 4, the second term corresponds to aerosol emissions due to volcanic eruptions of VEI greater than 3, and k is the decay rate of aerosols. A/12 should correspond roughly to emissions of volcanic eruptions of VEI less than 4 for 1 month, and 1 - (1-k)^(1/12) should give the rate of decay of volcanic aerosols between consecutive months. I get A = 0.0036 and k = 0.4758. Using the fact that there is no volcanic activity of VEI 4 or higher in 2012 and that the optical aerosol thickness in September 2012 is 0.0035, I get that the optical aerosol thickness in October, November and December should be approximately 0.0036, 0.0037 and 0.0038 respectively. LOD: I get length of day data from http://www.iers.org/IERS/EN/Science/EarthRotation/LODsince1623.html?nn=12932. As before, I lag everything by 6 years and detrend the data by 1.7 ms per century. SOI: http://www.bom.gov.au/climate/current/soihtm1.shtml PDO: http://www.ncdc.noaa.gov/teleconnections/pdo/data.csv I am a bit concerned that there might be a long term trend in the PDO due to global warming, but it is hard to tell. NAO: https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-station-based

I tried to use an a priori justified initial guess to see if it would lead to the gauss-newton estimation converging to something more reasonable: - If the fast feedback response with decay time 0.5 years of a doubling of CO2 concentrations causes temperatures to rise 1.15 C, then this would suggest that γ1 is 1.15/5.35/ln(2) = 0.3101. - If the remaining 4 exponential responses result in a climate sensitivity of approximately 3 C (and they are all of equal strength), then this would suggest that γ2 = γ3 = γ4 = γ5 = (3-1.15)/4/5.35/ln(2) = 0.1247. - The solar irradiance data is given in W/m^2. However, the Earth has an albedo of ~0.35, and the surface area of a sphere divided by the circle of the same radius is 4. Therefore, the coefficient solar in my model should be approximately 1/4*0.65 = 0.1625. - According to http://data.giss.nasa.gov/modelforce/strataer/, one can convert the aerosol optical depth to a forcing in W/m^2 by multiplying it by a factor of -23. Therefore, the volcano coefficient should be approximately -23. - I performed a simple exponential decay to equilibrium model and got coefficients for each of the internal variability indices. This coefficients should be roughly the same in the more complicated model. So I tried the above 'guess' with dG = 0. It still didn't converge to anything realistic either. - - - - - A second thing I thought I would try is to integrate my model so that temperature is the dependant variable rather than change in temperature. By integrating, my model becomes: The difference in temperature between January 1959 and month f + 1 is: Sum(s = 1 to 5; γs*((Sum i = 1 to f; dGHGi*ρs(i,i) + ... + dGHG1*ρs(i,1)) + Solar*(Sum i = 1 to f; dSi*ρs(i,i) + ... + dS1*ρs(i,1)) + Volcano*(Sum i = 1 to f; dVi*ρs(i,i) + ... + dV1*ρs(i,1)) + dGs*(Sum i = 1 to f; ρs(i,0)))) + β1*(LODf+1 – LOD1) + β2*(AAMf+1 – AAM1) + β3*(SOIf+1 – SOI1) + β4*(PDOf+1 – PDO1) + β5*(NAOf+1 – NAO1)) where month 1 corresponds to January 1959 Estimating this model still yields nonsense. The gauss-newton method stopped converging after only 1 iteration. - - - - - So the next thing I thought I would do is try to estimate the above model by using a sequence of linear regressions (alternating between keeping Solar, Volcano and dGs fixed and keeping the γs fixed). This gave slightly better results, but still somewhat meaningless. The next thing I tried was perform the sequence of linear regressions, but keep going (rather than stopping once the error stops decreasing). What I find is that as the number of iterations approach infinity, the γs approach zero and dG gets large. Obviously zero climate sensitivity is unrealistic. It seems like my attempts at trying to estimate the model are not converging to a reasonable solution. It might not be possible to estimate the model correctly with the 1959-2012 data set and it seems like the dG terms tend to soak up more of the variation. One option is that I try to use a different data set (say annual data from 1876 to 2011) and hope for the best; though data quality would decrease and some indices like AAM don't go that far back. A second option is to start from a period of time where it is safe for me to assume that dG = 0 (such as say 1750); unfortunately most data sets don't go back far enough to do this. A third option is the initially assume dG = 0, estimate the model for 1959-2012, and then use the estimated parameters to computer dG in 1959 by using annual data prior to 1959 (and repeating until convergence is met).

Also, I thought of another possible explanation for my nonsense results. The change in temperature is my dependent variable. Even if I do model the change in temperature very well, once I integrate to get temperature, the resulting temperature will be subject to 'random walks' away from the true temperature. That is the case if the assumptions about the error terms in the regression model are correct (independent, no heteroskedasticity, etc.). In reality, the behaviour of the error varies over time (example: the quality of the data in the 60's is not as good as the quality of the data in the 90's). When I take the residual of my results and I integrate, it is clear that what I get is subject to this random walk behaviour, which can explain my nonsense results. Maybe I should just integrate my current regression equation to get a new regression equation that has temperature as the dependent variable rather than change in temperature? Another thing I noticed is that rather than doing the Gauss-Newton method, I could also estimate the non-linear model by performing a sequence of linear regressions where I alternate between holding the impulse response function and the relative strength of different types of radiative forcings constant. This might give more robust results that is less dependent on an initial guess.

Thank you for your input Bonam. A few reasons. Firstly, I want enough data to be able to estimate the model and good data starts in 1959. The Mauna Loa CO2 data starts around this time (if I want CO2 data before this, I would have to use Law Dome ice core data which is a lot more questionable). The atmospheric angular momentum, southern oscillation index, solar magnetic flux and north atlantic multidecadal oscillation index data sets start in the 50's. In some cases there are reconstructions, but those reconstructions are more questionable.Also, before the 50's, temperature data becomes a lot more questionable. The temperature data during the world wars isn't very good and the global temperature data from this time period has been adjusted many times. Secondly, if I want to get a good estimate of the magnitude of response with a half year decay time then I probably want data with a higher temporal resolution that half a year. This is especially true if there are events like Pinotubo. Lastly, most of the data sets are given by month, so monthly data makes sense. Of course one has to take into account the unequal number of days, but that isn't impossible. Though I might want to try the model using annual data from say 1876-2012. One thing I was thinking of doing is to use multiple time series data sets (say monthly from 1959-2012, annual from 1876-1958, annual from 1650-1876, and the entire Pleistocene data over the past 400,000 years) and estimate the parameters for each time period. The idea is that the later time periods would be better at estimating the responses with short decay times, where as the earlier time periods would be better at estimating the responses with longer decay times. Provided the data sets are independent, then one could combine the results over all time periods to get a decent estimate of the impulse response function. Of course, it is better to resolve the current issues I am having before complicating things. Sorry, that was a typo. I wouldn't say it has no physical justification. Certainly the 0.5 year decay time has a physical justification. And since many responses to a small change in forcing should have a roughly exponential response (this would follow if the rate of change of the response is proportional to the magnitude of displacement from equilibrium) it probably is reasonable to expect that the impulse response function be a sum of exponentials. The impulse response function should be Integral(τ = 0 to infinity; f(τ)*(1-exp(-t/τ))dτ), where f(τ) is some unknown function of the decay time. The problem is that f(τ) is unknown. With respect to using decay time constants that are based on real times associated with known physical phenomena, I would like to do that. I tried to get a characteristic decay time of the Earth's heat sink towards equilibrium for example. But even in that, it proved difficult. In the case of the decay time of the Earth's heat sink, different heat sinks will have different decay times. A shallow ocean will have a faster response than a deeper ocean. And the surface of the ocean may have a faster response than the bottom of the ocean. As a result, rather than a single decay time, one gets a continuum of decay times with a large spread. So f(τ), rather than consisting only of a finite dirac deltas, is probably continuous. With respect to using a full orthonormal set to describe any function in that space, there are an infinite number of them (1-exp(-t/τ)), where τ is any real number between 0 and infinity. So you cannot estimate them all. And no clear first n orthonormal terms exist. So the next best thing that you can do is have a finite set that is distributed over all of the relevant values of τ (say from 0.5 years to 128 years). The other thing is that you probably want a higher density of these exponential decay functions where you think there is a higher density of f(τ). In our case, the smaller values of τ are generally going to be more relevant; for example, if one doubles CO2, then one should expect that the 0.5 year decay time term results in an increase of global temperatures of about 1.15K, then the next fastest term might be the water vapour feedback, which results in an increase of global temperatures of say 1K, then one has the snow-albedo feedback which might have a response time of about a decade, then one has the ocean heat sink response that has a response time of about 6 decades and brings to total change in temperatures to the ECS of about 3K, then one has the millennia long responses in vegetation and ice sheets which bring the total change in temperature up to the Earth system sensitivity of about 4.2K, etc. This justifies having a greater representation of small values of τ compared to larger values of τ (so having the finite choice of τ geometrically spaced over the relevant values of τ makes might be okay). So the Van Hateren approach is just a way to numerically determine the impulse response function from the data if you don't have much more a priori information about f(τ). Although a factor of 4 between consecutive τ might be too large. If I added terms with decay times of 1 year, 4 years, 16 years and 64 years, that would reduce the difference between consecutive τ to a factor of 2, which might be better. Given that my remaining degrees of freedom is over 600, reducing degrees of freedom by 4 would not be a big deal. Of course, this would make the individual estimates of the magnitudes of each exponential decay function less significant. But since I am interested primarily in the resulting impulse response function (and all the exponential response functions are strongly correlated with each other), this shouldn't result in any significant increase in uncertainty about the estimated impulse response function. I am not trying to use the exponential response functions to describe cyclical variability. I am directly using indices such as LOD, AAM, SOI, PDO and NAO to account for cyclical variability. The exponential response functions are supposed to describe how the climate moves towards a new equilibrium given a change in radiative forcing (be it more greenhouse gases, more solar irradiance, or more volcanic aerosols). I agree that sinusoidal functions would be a good way to represent cyclical climate variability. However, if I already have data on factors associated with climate variability (LOD, AAM, SOI, PDO and NAO), then why not use them directly? As an aside, I think that events like El-Nino are too chaotic to describe with sinusoidal functions, and that one might need to use Mathieu functions to represent them (if I were going to do that). Hyperbolics, would describe runaway phenomena, but the problem is that runaway phenomena are completely unrealistic for the Earth's climate response. Runaway global warming for Earth is basically impossible until the Earth reaches a temperature of 647 K. http://climatephys.org/2012/07/31/the-water-vapor-feedback-and-runaway-greenhouse/