<jupyter_start><jupyter_text>Problem Set 1<jupyter_code># imports
import numpy as np
from scipy.ndimage import uniform_filter1d
from matplotlib import pyplot as plt
import pandas
from datetime import datetime, timedelta

from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.api as sm
def set_fontsize(ax,fsz):
    """
    Set the fontsize throughout an Axis
    
    Args:
        ax (Matplotlib Axis): 
        fsz (float): Font size

    Returns:

    """
    for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
                 ax.get_xticklabels() + ax.get_yticklabels()):
        item.set_fontsize(fsz)<jupyter_output><empty_output><jupyter_text>Question 1  Load data<jupyter_code>data = pandas.read_table('../Data/greenland_mass_200204_202201.txt', 
                         delim_whitespace=True,
                         skiprows=31,
                       names=['decyear', 'mass', 'emass'])
data.head()<jupyter_output><empty_output><jupyter_text>Datetime<jupyter_code>dates = []

for start in data.decyear:
    #start = 2002.29
    year = int(start)
    rem = start - year

    base = datetime(year, 1, 1)
    result = base + timedelta(seconds=(base.replace(year=base.year + 1) - base).total_seconds() * rem)
    dates.append(result)
dates[0:5]
data['date'] = dates
data.set_index('date', inplace=True)
data.head()<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
ax = data.mass.plot(marker='o', ls='', ms=2)
# Label
ax.set_ylabel('Mass')
set_fontsize(ax,17)
plt.savefig('Problem1_timeseries.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Problem 2 -- Decompose<jupyter_code>result = seasonal_decompose(data.mass, model='additive', period=12)<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>fig = result.plot()
fig.set_size_inches((14, 11))
# Tight layout to realign things
fig.tight_layout()
plt.savefig('Problem2_decompose.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Problem 3 -- Description The trend is relatively linear, it is not grossly exponential. There is a strong seasonal component of approximaely 150 gt variation, peak-to-peak over the course of a year.  There is also an unusual 'hitch' during the first few months of the year.  There is an apparent seasonal variation to the residuals too, i.e. these two components "interact" The variance in the remainder increases during the first ~10 years and then appears relatively constant afterwards. ----  Problem 4 -- Berkeley monthly  Load<jupyter_code>globalT = pandas.read_table('../Data/LOT_berkeley_1850_2021_monthly.txt',
                         delim_whitespace=True,
                         skiprows=85, 
                         names = ['Year', 'Month',  'AnomalyM', 'Unc.M',   'AnomalyA', 'Unc.A',   
                                  'Anomaly5', 'Unc.5',   'Anomaly10',  'Unc.10',   'Anomaly20', 'Unc.20']
                           )
globalT.head()<jupyter_output><empty_output><jupyter_text>Dates<jupyter_code>dates = []
for index, row in globalT.iterrows():
    dates.append(f'{int(row.Year)}-{int(row.Month)}')
dates = pandas.to_datetime(dates)
globalT['date'] = dates
globalT.set_index('date', inplace=True)
globalT.head()<jupyter_output><empty_output><jupyter_text>Time series<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
ax = globalT.AnomalyM.plot(marker='o', ls='', ms=2)
# Label
ax.set_ylabel('Monthly Anomaly [K]')
set_fontsize(ax,17)
#plt.savefig('Problem1_timeseries.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>10 year filter<jupyter_code>filt_data = uniform_filter1d(globalT.AnomalyM, size=12*10)
globalT['filter10'] = filt_data
globalT.filter10<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
# Time series
ax = globalT.AnomalyM.plot(marker='o', ls='', ms=2)

# 10 year
globalT.filter10.plot(ax=ax, label='10 Year Filter')

# Label
ax.set_ylabel('Monthly Anomaly [K]')
set_fontsize(ax,17)
ax.legend(fontsize=17)
plt.savefig('Problem4_10year.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Problem 5 -- Other filters  Filter<jupyter_code>for nyear in [2,5,20]:
    filt_data = uniform_filter1d(globalT.AnomalyM, size=12*nyear)
    globalT[f'filter{nyear}'] = filt_data
globalT.head()<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
# Time series
ax = globalT.AnomalyM.plot(marker='o', ls='', ms=2)

# Filters
for nyear in [2,5,10,20]:
    getattr(globalT, f'filter{nyear}').plot(ax=ax, label=f'{nyear} Year Filter')

# Label
ax.set_ylabel('Monthly Anomaly [K]')
set_fontsize(ax,17)
ax.legend(fontsize=17)
plt.savefig('Problem5_multiyear.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>I would use 10 or 20 years to capture the overall trend of the data.   Shorter filters are too sensitive to variations that are on top of the overall trend. I think I favor 10 years but without a quantitative reason to prefer it. ----  Problem 7 -- Linear trend Arctice sea-ice  Load<jupyter_code>arctic_ice = pandas.read_csv('../Data/Sept_Arctic_extent_1979-2021.csv')
arctic_ice.head()<jupyter_output><empty_output><jupyter_text>Setup<jupyter_code>x = arctic_ice.year.values
X = sm.add_constant(x)
Y = arctic_ice.extent.values<jupyter_output><empty_output><jupyter_text>Fit<jupyter_code>arctic_model = sm.OLS(Y, X).fit()
predictions = arctic_model.predict(X) 
arctic_model.summary()<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
ax = plt.gca()
ax.plot(x, Y, "o", label="data")
#ax.plot(x, y_true, "b-", label="True")
ax.plot(x, arctic_model.fittedvalues, "r--.", label="OLS")
#ax.plot(x, iv_u, "r--")
#ax.plot(x, iv_l, "r--")
ax.legend(loc="best", fontsize=17)
ax.set_xlabel('Year', fontsize=17)
ax.set_ylabel('Extent (million sq km)', fontsize=17)
set_fontsize(ax, 17.)
#
plt.savefig('Problem7_arctic_linear.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Problem 8 -- What is the trend?<jupyter_code>arctic_ice.year.values[-1] - arctic_ice.year.values[0]
42 * arctic_model.params[1]
slope = arctic_model.params[1]
slope # million sq km / year<jupyter_output><empty_output><jupyter_text>Problem 9 -- Residuals<jupyter_code>arctic_resid = arctic_model.fittedvalues-Y<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>plt.clf()
fig = plt.figure(figsize=(12,8))
ax = plt.gca()
#ax.plot(x, Y, "o", label="data")
#ax.plot(x, y_true, "b-", label="True")
ax.plot(x, arctic_resid, "o", label="Residuals")
#ax.plot(x, iv_u, "r--")
#ax.plot(x, iv_l, "r--")
ax.legend(loc="best", fontsize=17)
ax.set_xlabel('Year', fontsize=17)
ax.set_ylabel('Extent (million sq km)', fontsize=17)
set_fontsize(ax, 17.)
plt.show()<jupyter_output><empty_output><jupyter_text>Normal?<jupyter_code>rms = np.std(arctic_resid)
rms
bin_size = 0.25
area = arctic_resid.size * bin_size<jupyter_output><empty_output><jupyter_text>Gaussian<jupyter_code>xval = np.linspace(-1.5,1.5,100)
yval = area * np.exp(-xval**2 / 2 / rms**2) / np.sqrt(2*np.pi*rms)
plt.clf()
fig = plt.figure(figsize=(10,6))
ax = plt.gca()
ax.hist(arctic_resid, bins=-1.5 + bin_size*np.arange(12))
# Overplot a Gaussian
ax.plot(xval, yval, color='r')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15234588106928532, 'fraction_numerical': 0.02905073649754501, 'mean_word_length': 5.289022298456261, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 89, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '543304', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3029, 'n_tokens_neox': 2825, 'n_words': 501, 'path': 'OCEA-267/ProblemSets/ProblemSet1.ipynb', 'repo_name': 'profxj/ocea200', 'size': '7332'}