As greenhouse gas we consider only CO$_2$. But the increase of CO$_2$ in the atmosphere should correlate strongly with the emission of other greenhouse gases.
# data from http://www.esrl.noaa.gov/gmd/ccgg/trends/
# since 1959 direct measurements!
co2_concentration = pd.read_csv("ghg-concentrations_fig-1_.csv")
# set year as index
co2_concentration = co2_concentration.set_index("Year")
co2_concentration.tail(3)
EPICA Dome C and Vostok Station, Antarctica | Law Dome, Antarctica (75-year smoothed) | Siple Station, Antarctica | Mauna Loa, Hawaii | Barrow, Alaska | Cape Matatula, American Samoa | South Pole, Antarctica | Cape Grim, Australia | Lampedusa Island, Italy | Shetland Islands, Scotland | |
---|---|---|---|---|---|---|---|---|---|---|
Year | ||||||||||
2012 | NaN | NaN | NaN | 393.82 | 394.940000 | 391.568333 | 390.060000 | NaN | NaN | NaN |
2013 | NaN | NaN | NaN | 396.48 | 397.928333 | 394.257500 | 392.779167 | NaN | NaN | NaN |
2014 | NaN | NaN | NaN | 398.55 | NaN | NaN | NaN | NaN | NaN | NaN |
Use Pandas for computing the mean of each row. NaN values are ignored:
# set new mean series
co2_concentration['mean'] = co2_concentration.mean(axis=1)
co2_concentration.tail(3)
EPICA Dome C and Vostok Station, Antarctica | Law Dome, Antarctica (75-year smoothed) | Siple Station, Antarctica | Mauna Loa, Hawaii | Barrow, Alaska | Cape Matatula, American Samoa | South Pole, Antarctica | Cape Grim, Australia | Lampedusa Island, Italy | Shetland Islands, Scotland | mean | |
---|---|---|---|---|---|---|---|---|---|---|---|
Year | |||||||||||
2012 | NaN | NaN | NaN | 393.82 | 394.940000 | 391.568333 | 390.060000 | NaN | NaN | NaN | 392.597083 |
2013 | NaN | NaN | NaN | 396.48 | 397.928333 | 394.257500 | 392.779167 | NaN | NaN | NaN | 395.361250 |
2014 | NaN | NaN | NaN | 398.55 | NaN | NaN | NaN | NaN | NaN | NaN | 398.550000 |
plt.plot(co2_concentration.index, co2_concentration['mean'])
plt.xlim((1400, 2020))
plt.xlabel("Year")
plt.ylabel("CO$_2$ concentration / ppm")
<matplotlib.text.Text at 0x10b3b9e50>
# convert c02 concentation to numpy
co2_year = co2_concentration.index.get_values()
co2_mean = co2_concentration['mean'].get_values()
lean_df = pd.read_csv("./lean2000_irradiance_only_data.csv", sep=',')
#lean_df
# solar irridiance
year_mask = (lean_df["YEAR"]- 0.5 >= df_air.Year.min())
scy = "11yrCYCLE+BKGRND"
#scy = "11yrCYCLE"
plt.plot(lean_df["YEAR"][year_mask],lean_df[scy][year_mask], "b*")
smo_ = np.ndarray(len(lean_df.YEAR))
import pandas.stats.moments
smoothed = pandas.stats.moments.ewma(lean_df[scy], halflife=3.)
plt.plot(lean_df["YEAR"][year_mask],smoothed[year_mask], "g-")
l = len(smoothed[year_mask])
smo_[:l] = smoothed[year_mask]
# according to a graph in http://www.climate4you.com/Sun.htm
# the irridience goes sligtly down sice 2000, to be more conservative we leave it constant
# from 2000-2014
smo_[l:] = smoothed[year_mask].max()
plt.xlabel("Year")
plt.ylabel("Solar irradiance / ?")
/Users/christian/.virtualenvs/alternative_1/lib/python2.7/site-packages/ipykernel/__main__.py:10: FutureWarning: pd.ewm_mean is deprecated for Series and will be removed in a future version, replace with Series.ewm(halflife=3.0,ignore_na=False,min_periods=0,adjust=True).mean()
<matplotlib.text.Text at 0x10b59ee90>
#mean = (df_air["J-D"]*0.01).mean()
temp_dev_air = -0.2 + ((df_air["J-D"]-df_air["J-D"].min())*0.01).get_values()
#temp_dev_lo = ((df_land_and_ocean["J-D"] - (df_land_and_ocean["J-D"]).min())*0.01).get_values()
plot_tempe()
plt.plot(co2_original, alo_temp, 'bo')
plt.xlabel("CO$_2$ / ppm")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")
<matplotlib.text.Text at 0x10b5f1dd0>
l = (year_mask).sum()
plt.plot(smoothed[year_mask], alo_temp[:l], 'bo')
plt.xlabel("solar irridiance / ?")
plt.ylabel("$\Delta$ temperature / ${}^\circ$C")
<matplotlib.text.Text at 0x10ba22ed0>
import pymc3
model = pymc3.Model()
with model:
co_2_weight = pymc3.Normal("co_2_weight", mu=0, sd=1.)
irridiation_weight = pymc3.Normal("irridiaton_weight", mu=0, sd=1.)
offset = pymc3.Normal("offset", mu=0, sd=1.)
i_year = year - year[0]
i_co2_year = year - co2_year[0] # + co_2_delay
temp_mu = offset + irridiation_weight * smo[i_year] + co_2_weight * co2_
temp_sigma = pymc3.HalfNormal('sigma', sd=0.5)
temp_dev_ = pymc3.Normal("temp_dev_", mu=temp_mu, sd=temp_sigma, observed = alo_temp)
temp_pred = pymc3.Normal("temp_pred", mu=temp_mu, sd=temp_sigma, shape = alo_temp.shape )
Couldn't import dot_parser, loading of dot files will not be possible. Applied log-transform to sigma and added transformed sigma_log to model.
with model:
# obtain starting values via MAP
start = pymc3.find_MAP()#fmin=optimize.fmin_powell)
#step = pymc3.Slice()
tr = pymc3.sample(5000, start=start)
Assigned NUTS to co_2_weight Assigned NUTS to irridiaton_weight Assigned NUTS to offset Assigned NUTS to sigma_log Assigned NUTS to temp_pred [-----------------100%-----------------] 5000 of 5000 complete in 7.2 sec
plt.figure(figsize=(7, 7))
pymc3.traceplot(tr)
plt.tight_layout();
<matplotlib.figure.Figure at 0x1143897d0>
Compare the influence of irridiation and CO$_2$(Normalized Data):
irridiation_weight_.mean()
0.045895907084405878
co_2_weight_.mean()
0.33630017388761918
plt.figure(figsize=(12,10))
temp_pred_ = tr.get_values(temp_pred)[4000:]
for t in temp_pred_:
plt.plot(year, t, 'b.', alpha=0.02)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')
[<matplotlib.lines.Line2D at 0x113cc7b10>]
plt.figure(figsize=(12,10))
for t in temp_.T:
plt.plot(year, t, 'b-', alpha=0.04)
plt.xlabel("Year")
plt.ylabel("$\Delta T$ / ${}^\circ$C")
plt.plot(year, alo_temp, 'k.')
[<matplotlib.lines.Line2D at 0x113cc7b90>]
import sklearn
import sklearn.linear_model
lr = sklearn.linear_model.LinearRegression()
lr.fit(last_years[:,None], co2_last_years)
future_co2 = lr.predict(future[:,None])
/Users/christian/.virtualenvs/alternative_1/lib/python2.7/site-packages/scipy/linalg/basic.py:884: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver. warnings.warn(mesg, RuntimeWarning)
Assumption: The CO$_2$ level in atmosphere will (linear) increase as in the last years (from 1980 to now).
with pymc3.Model() as model_co2:
intercept = pymc3.Normal("intercept", mu=0, sd=100.)
l_weight = pymc3.Normal("l_weight", mu=0, sd=10.)
co2_mu = intercept + l_weight * last_years_shared
co2_sigma = pymc3.HalfNormal('co2_sigma', sd=0.2)
co2_pred = pymc3.Normal("co2_pred", mu=co2_mu, sd=co2_sigma, observed = co2_last_years)
co2_sample = pymc3.Normal("co2_sample", mu=co2_mu, sd=co2_sigma, shape = co2_last_years.shape)
Applied log-transform to co2_sigma and added transformed co2_sigma_log to model.
with model_co2:
#step = pymc3.Slice()
trace_co2 = pymc3.sample(10000)
Assigned NUTS to intercept Assigned NUTS to l_weight Assigned NUTS to co2_sigma_log Assigned NUTS to co2_sample [-----------------100%-----------------] 10000 of 10000 complete in 11.2 sec
plt.figure(figsize=(7, 7))
pymc3.traceplot(trace_co2)
plt.tight_layout();
<matplotlib.figure.Figure at 0x12d7c3610>
nb_samples=500
last_years_shared.set_value(future)
# Simply running PPC will use the updated values and do prediction
ppc = pymc3.sample_ppc(trace_co2[5000:], model=model_co2, samples=nb_samples)
co2_future = ppc["co2_pred"]
co2_future.shape
(500, 40)
plt.figure(figsize=(14,10))
for t in temp_[-len(temp_future_.T):].T:
plt.plot(year, t, 'b.', alpha=0.02)
for t in temp_future_.T:
plt.plot(future, t, 'r.', alpha=0.02)
plt.plot(year, alo_temp, 'k.')
plt.ylabel("$\delta T$ / C")
plt.xlabel("Year")
<matplotlib.text.Text at 0x145b1f250>