app weather innovation

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%22marimo%3E%3D0.23.9%22%2C%20%22matplotlib%3D%3D3.10.8%22%2C%20%22numpy%3D%3D2.4.3%22%2C%20%22pandas%3D%3D3.0.2%22%2C%20%22scipy%3D%3D1.17.1%22%5D%0A%23%20%0A%23%20%5Btool.marimo.export%5D%0A%23%20lock_kind%20%3D%20%22resolved%22%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.1%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Berkeley%20Earth%20Temperature%20Anomalies%3A%20Innovation%20Filtering%20via%20Spectral%20Factorization%0A%0A%20%20%20%20%23%23%20Overview%0A%0A%20%20%20%20This%20notebook%20demonstrates%20how%20to%20extract%20%22innovations%22%20(white%20noise%20shocks)%20from%20temperature%0A%20%20%20%20anomalies%20(daily%20or%20monthly)%20using%20spectral%20factorization%20and%20cepstral%20analysis.%20The%20innovations%20represent%20the%0A%20%20%20%20unpredictable%20component%20of%20weather%20variability%2C%20while%20the%20autocorrelation%20structure%20captures%0A%20%20%20%20persistent%20patterns%20(e.g.%2C%20multi-day%20or%20multi-month%20weather%20systems).%0A%0A%20%20%20%20**Key%20concepts%3A**%0A%0A%20%20%20%20-%20Temperature%20anomalies%20from%20Berkeley%20Earth%20represent%20deviations%20from%20climatological%20norms%0A%20%20%20%20-%20Autocorrelation%20and%20Power%20Spectral%20Density%20(PSD)%20characterize%20the%20temporal%20structure%0A%20%20%20%20-%20Cepstral%20spectral%20factorization%20yields%20a%20minimum-phase%20shaping%20filter%20H(%CF%89)%20such%20that%0A%20%20%20%20%20%20S_x(%CF%89)%20%E2%89%88%20%CF%83%C2%B2%7CH(%CF%89)%7C%C2%B2%0A%20%20%20%20-%20The%20innovation%20(whitening)%20filter%20G(%CF%89)%20%3D%201%2FH(%CF%89)%20extracts%20white%20innovations%20from%20the%20anomalies%0A%20%20%20%20-%20These%20innovations%20can%20be%20used%20for%20one-step-ahead%20prediction%20via%20AR%20modeling%0A%0A%20%20%20%20**Narrative%3A**%20Weather%20anomalies%20exhibit%20persistence%20(autocorrelation)%2C%20suggesting%20they%20are%20the%0A%20%20%20%20output%20of%20a%20linear%20system%20driven%20by%20white%20%22weather%20shocks.%22%20The%20innovation%20filter%20inverts%20this%0A%20%20%20%20system%20to%20recover%20those%20shocks%2C%20which%20are%20approximately%20white%20and%20uncorrelated.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20scipy%20import%20signal%0A%20%20%20%20from%20scipy%20import%20fft%0A%20%20%20%20import%20warnings%0A%20%20%20%20warnings.filterwarnings('ignore')%0A%0A%20%20%20%20%23%20Optional%3A%20statsmodels%20for%20AR%20model%20fitting%20(for%20comparison)%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20import%20statsmodels.api%20as%20sm%0A%20%20%20%20%20%20%20%20from%20statsmodels.tsa.stattools%20import%20acf%20as%20sm_acf%0A%20%20%20%20%20%20%20%20from%20statsmodels.tsa.ar_model%20import%20AutoReg%0A%20%20%20%20%20%20%20%20HAS_STATSMODELS%20%3D%20True%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20HAS_STATSMODELS%20%3D%20False%0A%20%20%20%20%20%20%20%20print(%22Note%3A%20statsmodels%20not%20available.%20AR%20fitting%20will%20use%20Yule-Walker%20equations%20only.%22)%0A%0A%20%20%20%20%23%20Configuration%0A%20%20%20%20np.random.seed(42)%20%20%23%20For%20reproducibility%0A%20%20%20%20plt.style.use('default')%0A%20%20%20%20plt.rcParams%5B'figure.figsize'%5D%20%3D%20(12%2C%206)%0A%20%20%20%20plt.rcParams%5B'font.size'%5D%20%3D%2010%0A%20%20%20%20return%20AutoReg%2C%20HAS_STATSMODELS%2C%20fft%2C%20np%2C%20pd%2C%20plt%2C%20signal%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt)%3A%0A%20%20%20%20def%20load_berkeley_temperature(path%3A%20str)%20-%3E%20pd.Series%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Load%20Berkeley%20Earth%20temperature%20CSV%20(daily%20or%20monthly)%20and%20return%20as%20a%20pandas%20Series%20with%20DatetimeIndex.%0A%0A%20%20%20%20%20%20%20%20Handles%20both%3A%0A%20%20%20%20%20%20%20%20-%20Standard%20format%3A%20CSV%20with%20Date%20and%20Temperature%20columns%20(typically%20daily)%0A%20%20%20%20%20%20%20%20-%20Monthly%20format%3A%20Berkeley%20Earth%20format%20with%20comment%20lines%20and%20Year%2FMonth%20columns%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20-----------%0A%20%20%20%20%20%20%20%20path%20%3A%20str%0A%20%20%20%20%20%20%20%20%20%20%20%20Path%20to%20the%20CSV%20file%0A%20%20%20%20%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20--------%0A%20%20%20%20%20%20%20%20pd.Series%0A%20%20%20%20%20%20%20%20%20%20%20%20Temperature%20series%20with%20DatetimeIndex%2C%20sorted%20by%20date%2C%20NaNs%20removed%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%20%20%23%20Try%20reading%20as%20space-separated%20first%20(Berkeley%20Earth%20monthly%20format)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20df%20%3D%20pd.read_csv(path%2C%20comment%3D'%25'%2C%20sep%3Dr'%5Cs%2B'%2C%20header%3DNone%2C%20skipinitialspace%3DTrue%2C%20na_values%3D%5B'NaN'%2C%20'nan'%5D)%20%20%23%20This%20format%20has%20comment%20lines%20starting%20with%20%25%20and%20space-separated%20data%20with%20no%20header%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20df.shape%5B1%5D%20%3E%3D%202%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20df.columns%20%3D%20%5B'Year'%2C%20'Month'%2C%20'Monthly_Anomaly'%2C%20'Monthly_Unc'%5D%20%2B%20%5Bf'Col_%7Bi%7D'%20for%20i%20in%20range(4%2C%20df.shape%5B1%5D)%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20df%20%3D%20df.dropna(subset%3D%5B'Year'%2C%20'Month'%5D%2C%20how%3D'all')%0A%20%20%20%20%20%20%20%20%20%20%20%20except%3A%20%20%23%20Check%20if%20this%20looks%20like%20Berkeley%20Earth%20monthly%20format%20(Year%2C%20Month%20columns)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20df%20%3D%20pd.read_csv(path%2C%20comment%3D'%25')%0A%20%20%20%20%20%20%20%20except%20FileNotFoundError%3A%20%20%23%20Assume%20columns%3A%20Year%2C%20Month%2C%20Monthly%20Anomaly%2C%20Monthly%20Uncertainty%2C%20...%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%20FileNotFoundError(f%22File%20not%20found%3A%20%7Bpath%7D%5Cn%5CnPlease%20download%20a%20Berkeley%20Earth%20temperature%20CSV%20file%20from%3A%5Cnhttps%3A%2F%2Fberkeleyearth.org%2Fdata%2F%5CnSave%20it%20as%20'data%2Fberkeley_city_daily.csv'%20or%20'data%2Fberkeley_city_monthly.csv'%22)%0A%20%20%20%20%20%20%20%20if%20'Year'%20in%20df.columns%20and%20'Month'%20in%20df.columns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5B'Year'%5D%20%3D%20pd.to_numeric(df%5B'Year'%5D%2C%20errors%3D'coerce')%20%20%23%20Remove%20any%20completely%20empty%20rows%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5B'Month'%5D%20%3D%20pd.to_numeric(df%5B'Month'%5D%2C%20errors%3D'coerce')%0A%20%20%20%20%20%20%20%20%20%20%20%20df%20%3D%20df.dropna(subset%3D%5B'Year'%2C%20'Month'%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5B'Year'%5D%20%3D%20df%5B'Year'%5D.astype(int)%20%20%23%20Fall%20back%20to%20standard%20CSV%20format%20(typically%20daily%20data)%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5B'Month'%5D%20%3D%20df%5B'Month'%5D.astype(int)%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5B'Date'%5D%20%3D%20pd.to_datetime(df%5B%5B'Year'%2C%20'Month'%5D%5D.assign(Day%3D1))%0A%20%20%20%20%20%20%20%20%20%20%20%20date_col%20%3D%20'Date'%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20'Monthly_Anomaly'%20in%20df.columns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20'Monthly_Anomaly'%0A%20%20%20%20%20%20%20%20%20%20%20%20elif%20df.shape%5B1%5D%20%3E%3D%203%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20df.columns%5B2%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print(f%22Using%20column%20'%7Btemp_col%7D'%20(3rd%20column)%20as%20monthly%20anomaly%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20numeric_cols%20%3D%20%5Bcol%20for%20col%20in%20df.columns%20if%20df%5Bcol%5D.dtype%20in%20%5Bnp.number%5D%20and%20col%20not%20in%20%5B'Year'%2C%20'Month'%5D%5D%20%20%23%20Check%20if%20this%20is%20Berkeley%20Earth%20monthly%20format%20(has%20Year%20and%20Month%20columns)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20len(numeric_cols)%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20numeric_cols%5B0%5D%20%20%23%20Monthly%20format%3A%20ensure%20Year%20and%20Month%20are%20numeric%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print(f%22Using%20column%20'%7Btemp_col%7D'%20as%20temperature%2Fanomaly%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError('Could%20not%20find%20a%20temperature%2Fanomaly%20column')%20%20%23%20Remove%20rows%20where%20Year%20or%20Month%20is%20NaN%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5Btemp_col%5D%20%3D%20pd.to_numeric(df%5Btemp_col%5D%2C%20errors%3D'coerce')%0A%20%20%20%20%20%20%20%20%20%20%20%20print(f'Loaded%20monthly%20data%3A%20%7Blen(df)%7D%20months')%20%20%23%20Convert%20to%20integers%0A%20%20%20%20%20%20%20%20%20%20%20%20print(f'Note%3A%20Data%20is%20already%20anomalies%20(relative%20to%201951-1980%20baseline)')%0A%20%20%20%20%20%20%20%20%20%20%20%20df.attrs%5B'is_anomaly'%5D%20%3D%20True%0A%20%20%20%20%20%20%20%20else%3A%20%20%23%20Create%20date%20from%20Year%20and%20Month%0A%20%20%20%20%20%20%20%20%20%20%20%20date_cols%20%3D%20%5B'Date'%2C%20'date'%2C%20'DATE'%2C%20'time'%2C%20'Time'%2C%20'TIME'%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20date_col%20%3D%20None%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20col%20in%20date_cols%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20col%20in%20df.columns%3A%20%20%23%20Use%20Monthly%20Anomaly%20column%20(3rd%20column%2C%20index%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20date_col%20%3D%20col%20%20%23%20Columns%20are%3A%20Year%2C%20Month%2C%20Monthly_Anomaly%2C%20Monthly_Unc%2C%20...%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20date_col%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20col%20in%20df.columns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20df%5Bcol%5D.dtype%20%3D%3D%20'object'%3A%20%20%23%20Use%203rd%20column%20(index%202)%20which%20should%20be%20the%20monthly%20anomaly%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pd.to_datetime(df%5Bcol%5D.iloc%5B0%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20date_col%20%3D%20col%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%20%20%23%20Find%20first%20numeric%20column%20after%20Year%20and%20Month%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20except%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pass%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20date_col%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError('Could%20not%20find%20a%20date%20column.%20Expected%20one%20of%3A%20'%20%2B%20'%2C%20'.join(date_cols))%0A%20%20%20%20%20%20%20%20%20%20%20%20temp_cols%20%3D%20%5B'Temperature'%2C%20'temperature'%2C%20'TEMP'%2C%20'temp'%2C%20'AverageTemperature'%2C%20'avg_temp'%2C%20'T'%2C%20't'%2C%20'Anomaly'%2C%20'anomaly'%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20None%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20col%20in%20temp_cols%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20col%20in%20df.columns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20col%20%20%23%20Convert%20temperature%20column%20to%20numeric%2C%20handling%20NaN%20strings%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20temp_col%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20numeric_cols%20%3D%20df.select_dtypes(include%3D%5Bnp.number%5D).columns%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20len(numeric_cols)%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20temp_col%20%3D%20numeric_cols%5B0%5D%20%20%23%20Mark%20that%20this%20data%20is%20already%20anomalies%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print(f%22Using%20numeric%20column%20'%7Btemp_col%7D'%20as%20temperature%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError('Could%20not%20find%20a%20temperature%20column.%20Expected%20one%20of%3A%20'%20%2B%20'%2C%20'.join(temp_cols))%0A%20%20%20%20%20%20%20%20%20%20%20%20df%5Bdate_col%5D%20%3D%20pd.to_datetime(df%5Bdate_col%5D)%20%20%23%20Standard%20format%20(typically%20daily)%3A%20try%20to%20detect%20date%20column%0A%20%20%20%20%20%20%20%20series%20%3D%20pd.Series(data%3Ddf%5Btemp_col%5D.values%2C%20index%3Ddf%5Bdate_col%5D.values%2C%20name%3D'temperature')%0A%20%20%20%20%20%20%20%20series%20%3D%20series.dropna().sort_index()%0A%20%20%20%20%20%20%20%20if%20hasattr(df%2C%20'attrs')%20and%20df.attrs.get('is_anomaly'%2C%20False)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20series.attrs%20%3D%20%7B'is_anomaly'%3A%20True%7D%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20series.attrs%20%3D%20%7B'is_anomaly'%3A%20False%7D%0A%20%20%20%20%20%20%20%20return%20series%0A%20%20%20%20data_paths%20%3D%20%5B'data%2Fberkeley_london_monthly.csv'%2C%20'data%2Fberkeley_city_daily.csv'%2C%20'data%2Fberkeley_city_monthly.csv'%5D%0A%20%20%20%20temp_series%20%3D%20None%20%20%23%20Try%20to%20find%20a%20column%20that%20looks%20like%20dates%0A%20%20%20%20for%20path%20in%20data_paths%3A%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20temp_series%20%3D%20load_berkeley_temperature(path)%0A%20%20%20%20%20%20%20%20%20%20%20%20print(f'Successfully%20loaded%20data%20from%3A%20%7Bpath%7D')%0A%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20%20%20%20%20except%20FileNotFoundError%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20if%20temp_series%20is%20None%3A%0A%20%20%20%20%20%20%20%20raise%20FileNotFoundError('Could%20not%20find%20temperature%20data%20file.%20Please%20ensure%20one%20of%20these%20exists%3A%5Cn-%20data%2Fberkeley_london_monthly.csv%5Cn-%20data%2Fberkeley_city_daily.csv%5Cn-%20data%2Fberkeley_city_monthly.csv')%0A%20%20%20%20_freq%20%3D%20pd.infer_freq(temp_series.index)%0A%20%20%20%20_median_days%20%3D%20temp_series.index.to_series().diff().dt.days.median()%0A%20%20%20%20_is_monthly%20%3D%20bool((_freq%20and%20('M'%20in%20_freq%20or%20'MS'%20in%20_freq))%20or%20((temp_series.index.day%20%3D%3D%201).all()%20and%20_median_days%20%3E%3D%2028))%0A%20%20%20%20if%20_is_monthly%3A%0A%20%20%20%20%20%20%20%20print(f'Loaded%20%7Blen(temp_series)%7D%20monthly%20temperature%20records')%20%20%23%20Try%20to%20detect%20temperature%20column%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20print(f'Loaded%20%7Blen(temp_series)%7D%20daily%20temperature%20records')%0A%20%20%20%20print(f'Date%20range%3A%20%7Btemp_series.index.min()%7D%20to%20%7Btemp_series.index.max()%7D')%0A%20%20%20%20plt.figure(figsize%3D(14%2C%205))%0A%20%20%20%20plt.plot(temp_series.index%2C%20temp_series.values%2C%20'b-'%2C%20alpha%3D0.6%2C%20linewidth%3D0.5)%0A%20%20%20%20plt.xlabel('Date')%0A%20%20%20%20plt.ylabel('Temperature%20(%C2%B0C)')%0A%20%20%20%20title%20%3D%20'Raw%20Monthly%20Temperature%20Data'%20if%20_is_monthly%20else%20'Raw%20Daily%20Temperature%20Data'%0A%20%20%20%20plt.title(title)%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%20%20%23%20Try%20numeric%20columns%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20%23%20Load%20the%20data%0A%20%20%20%20%23%20Try%20monthly%20first%2C%20then%20daily%0A%20%20%20%20%23%20Determine%20if%20data%20is%20monthly%20or%20daily%20based%20on%20frequency%0A%20%20%20%20%23%20Plot%20raw%20temperature%0A%20%20%20%20plt.show()%20%20%23%20Parse%20dates%20%20%23%20Create%20series%20%20%23%20Remove%20NaNs%20and%20sort%20%20%23%20Store%20metadata%20about%20whether%20data%20is%20already%20anomalies%0A%20%20%20%20return%20(temp_series%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20pd%2C%20plt%2C%20temp_series)%3A%0A%20%20%20%20from%20scipy.interpolate%20import%20interp1d%0A%0A%20%20%20%20dates%20%3D%20temp_series.index%0A%20%20%20%20temperatures%20%3D%20temp_series.values%0A%20%20%20%20N%20%3D%20len(temperatures)%0A%20%20%20%20_freq%20%3D%20pd.infer_freq(dates)%0A%20%20%20%20_median_days%20%3D%20dates.to_series().diff().dt.days.median()%0A%20%20%20%20_is_monthly%20%3D%20bool((_freq%20and%20('M'%20in%20_freq%20or%20'MS'%20in%20_freq))%20or%20((dates.day%20%3D%3D%201).all()%20and%20_median_days%20%3E%3D%2028))%0A%20%20%20%20is_already_anomaly%20%3D%20False%0A%20%20%20%20if%20_is_monthly%3A%0A%20%20%20%20%20%20%20%20mean_temp%20%3D%20np.mean(temperatures)%0A%20%20%20%20%20%20%20%20if%20abs(mean_temp)%20%3C%202.0%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20is_already_anomaly%20%3D%20True%0A%20%20%20%20%20%20%20%20%20%20%20%20print(f'Detected%20that%20data%20is%20already%20anomalies%20(mean%20%E2%89%88%20%7Bmean_temp%3A.4f%7D%20%C2%B0C)')%0A%20%20%20%20%20%20%20%20%20%20%20%20print('Skipping%20climatological%20cycle%20removal%20-%20only%20ensuring%20zero%20mean')%0A%20%20%20%20if%20_is_monthly%20and%20(not%20is_already_anomaly)%3A%0A%20%20%20%20%20%20%20%20month%20%3D%20dates.month.values%0A%20%20%20%20%20%20%20%20period_unique%20%3D%20np.arange(1%2C%2013)%0A%20%20%20%20%20%20%20%20clim_mean%20%3D%20np.zeros(12)%0A%20%20%20%20%20%20%20%20for%20m%20in%20range(1%2C%2013)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20mask%20%3D%20month%20%3D%3D%20m%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20np.any(mask)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%5Bm%20-%201%5D%20%3D%20np.mean(temperatures%5Bmask%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%5Bm%20-%201%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20valid_mask%20%3D%20~np.isnan(clim_mean)%0A%20%20%20%20%20%20%20%20if%20np.any(valid_mask)%20and%20np.any(~valid_mask)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20clim_interp%20%3D%20interp1d(period_unique%5Bvalid_mask%5D%2C%20clim_mean%5Bvalid_mask%5D%2C%20kind%3D'linear'%2C%20fill_value%3D'extrapolate'%2C%20bounds_error%3DFalse)%0A%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%20%3D%20clim_interp(period_unique)%0A%20%20%20%20%20%20%20%20clim_series%20%3D%20clim_mean%5Bmonth%20-%201%5D%0A%20%20%20%20%20%20%20%20period_label%20%3D%20'Month'%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20day_of_year%20%3D%20dates.dayofyear.values%0A%20%20%20%20%20%20%20%20period_unique%20%3D%20np.arange(1%2C%20367)%0A%20%20%20%20%20%20%20%20clim_mean%20%3D%20np.zeros(366)%0A%20%20%20%20%20%20%20%20for%20doy%20in%20range(1%2C%20367)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20mask%20%3D%20day_of_year%20%3D%3D%20doy%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20np.any(mask)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%5Bdoy%20-%201%5D%20%3D%20np.mean(temperatures%5Bmask%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%5Bdoy%20-%201%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20valid_mask%20%3D%20~np.isnan(clim_mean)%0A%20%20%20%20%20%20%20%20if%20np.any(valid_mask)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20clim_interp%20%3D%20interp1d(period_unique%5Bvalid_mask%5D%2C%20clim_mean%5Bvalid_mask%5D%2C%20kind%3D'linear'%2C%20fill_value%3D'extrapolate'%2C%20bounds_error%3DFalse)%0A%20%20%20%20%20%20%20%20%20%20%20%20clim_mean%20%3D%20clim_interp(period_unique)%0A%20%20%20%20%20%20%20%20clim_series%20%3D%20clim_mean%5Bday_of_year%20-%201%5D%0A%20%20%20%20%20%20%20%20period_label%20%3D%20'Day%20of%20Year'%0A%20%20%20%20if%20is_already_anomaly%3A%0A%20%20%20%20%20%20%20%20anomalies%20%3D%20temperatures%20-%20np.mean(temperatures)%0A%20%20%20%20%20%20%20%20print('Using%20data%20as-is%20(already%20anomalies)%2C%20only%20removed%20mean%20to%20ensure%20zero%20mean')%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20global_mean%20%3D%20np.mean(temperatures)%0A%20%20%20%20%20%20%20%20anomalies%20%3D%20temperatures%20-%20clim_series%0A%20%20%20%20%20%20%20%20anomalies%20%3D%20anomalies%20-%20np.mean(anomalies)%0A%20%20%20%20%20%20%20%20print('Removed%20climatological%20seasonal%20cycle%20and%20ensured%20zero%20mean')%0A%20%20%20%20anomaly_series%20%3D%20pd.Series(anomalies%2C%20index%3Ddates%2C%20name%3D'anomaly')%0A%20%20%20%20print(f'Anomaly%20statistics%3A')%0A%20%20%20%20print(f'%20%20Mean%3A%20%7Bnp.mean(anomalies)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20print(f'%20%20Std%3A%20%7Bnp.std(anomalies)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20print(f'%20%20Min%3A%20%7Bnp.min(anomalies)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20print(f'%20%20Max%3A%20%7Bnp.max(anomalies)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(14%2C%208)%2C%20sharex%3DTrue)%0A%20%20%20%20window_years%20%3D%203%0A%20%20%20%20if%20_is_monthly%3A%0A%20%20%20%20%20%20%20%20window_periods%20%3D%20window_years%20*%2012%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20window_periods%20%3D%20window_years%20*%20365%0A%20%20%20%20_start_idx%20%3D%20max(0%2C%20N%20-%20window_periods)%0A%20%20%20%20_axes%5B0%5D.plot(dates%5B_start_idx%3A%5D%2C%20temperatures%5B_start_idx%3A%5D%2C%20'b-'%2C%20alpha%3D0.7%2C%20linewidth%3D0.8%2C%20label%3D'Temperature'%20if%20is_already_anomaly%20else%20'Raw%20temperature')%0A%20%20%20%20if%20not%20is_already_anomaly%3A%0A%20%20%20%20%20%20%20%20_axes%5B0%5D.plot(dates%5B_start_idx%3A%5D%2C%20clim_series%5B_start_idx%3A%5D%2C%20'r--'%2C%20alpha%3D0.7%2C%20linewidth%3D1%2C%20label%3D'Climatological%20mean')%0A%20%20%20%20_axes%5B0%5D.set_ylabel('Temperature%20(%C2%B0C)')%0A%20%20%20%20title_label%20%3D%20'Temperature%20Anomalies'%20if%20is_already_anomaly%20else%20f'Raw%20Temperature%20and%20Climatological%20Mean%20(last%20%7Bwindow_years%7D%20years)'%0A%20%20%20%20_axes%5B0%5D.set_title(title_label)%0A%20%20%20%20_axes%5B0%5D.legend()%0A%20%20%20%20_axes%5B0%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20_axes%5B1%5D.plot(dates%5B_start_idx%3A%5D%2C%20anomalies%5B_start_idx%3A%5D%2C%20'g-'%2C%20alpha%3D0.7%2C%20linewidth%3D0.8)%0A%20%20%20%20_axes%5B1%5D.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20alpha%3D0.3)%0A%20%20%20%20_axes%5B1%5D.set_ylabel('Temperature%20Anomaly%20(%C2%B0C)')%0A%20%20%20%20_axes%5B1%5D.set_xlabel('Date')%0A%20%20%20%20_axes%5B1%5D.set_title('Temperature%20Anomalies')%0A%20%20%20%20_axes%5B1%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20if%20not%20is_already_anomaly%3A%0A%20%20%20%20%20%20%20%20plt.figure(figsize%3D(10%2C%205))%0A%20%20%20%20%20%20%20%20plt.plot(period_unique%2C%20clim_mean%2C%20'r-'%2C%20linewidth%3D2)%0A%20%20%20%20%20%20%20%20plt.xlabel(period_label)%0A%20%20%20%20%20%20%20%20plt.ylabel('Climatological%20Mean%20Temperature%20(%C2%B0C)')%0A%20%20%20%20%20%20%20%20cycle_title%20%3D%20'Monthly%20Climatological%20Cycle'%20if%20_is_monthly%20else%20'Annual%20Climatological%20Cycle'%0A%20%20%20%20%20%20%20%20plt.title(cycle_title)%0A%20%20%20%20%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20%20%20%20%20plt.tight_layout()%0A%20%20%20%20%20%20%20%20plt.show()%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20print('Skipping%20climatological%20cycle%20plot%20(data%20is%20already%20anomalies)')%0A%20%20%20%20return%20N%2C%20anomalies%2C%20dates%0A%0A%0A%40app.cell%0Adef%20_(anomalies%2C%20dates%2C%20np%2C%20pd%2C%20plt)%3A%0A%20%20%20%20def%20estimate_autocorrelation(x%3A%20np.ndarray%2C%20max_lag%3A%20int)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Estimate%20autocorrelation%20function%20using%20biased%20estimator.%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20-----------%0A%20%20%20%20%20%20%20%20x%20%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20Input%20signal%20(1D%20array)%0A%20%20%20%20%20%20%20%20max_lag%20%3A%20int%0A%20%20%20%20%20%20%20%20%20%20%20%20Maximum%20lag%20to%20compute%20(inclusive)%0A%20%20%20%20%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20--------%0A%20%20%20%20%20%20%20%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20Autocorrelation%20R%5B0..max_lag%5D%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20N%20%3D%20len(x)%0A%20%20%20%20%20%20%20%20R%20%3D%20np.zeros(max_lag%20%2B%201)%0A%0A%20%20%20%20%20%20%20%20for%20lag%20in%20range(max_lag%20%2B%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20lag%20%3C%20N%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20R%5Blag%5D%20%3D%20np.mean(x%5B%3AN-lag%5D%20*%20x%5Blag%3A%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20R%5Blag%5D%20%3D%200.0%0A%0A%20%20%20%20%20%20%20%20return%20R%0A%0A%0A%20%20%20%20%23%20Estimate%20autocorrelation%0A%20%20%20%20%23%20Adjust%20max_lag%20based%20on%20data%20frequency%0A%20%20%20%20freq_check%20%3D%20pd.infer_freq(dates)%0A%20%20%20%20median_days%20%3D%20dates.to_series().diff().dt.days.median()%0A%20%20%20%20is_monthly_data%20%3D%20bool((freq_check%20and%20('M'%20in%20freq_check%20or%20'MS'%20in%20freq_check))%20or%20((dates.day%20%3D%3D%201).all()%20and%20median_days%20%3E%3D%2028))%0A%20%20%20%20if%20is_monthly_data%3A%0A%20%20%20%20%20%20%20%20max_lag%20%3D%2024%20%20%23%20Two%20years%20of%20monthly%20lags%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20max_lag%20%3D%20365%20%20%23%20One%20year%20of%20lags%20(for%20daily%20data)%0A%20%20%20%20R_x%20%3D%20estimate_autocorrelation(anomalies%2C%20max_lag)%0A%20%20%20%20lags%20%3D%20np.arange(len(R_x))%0A%0A%20%20%20%20%23%20Normalize%20by%20R%5B0%5D%20to%20get%20correlation%20coefficient%0A%20%20%20%20R_x_normalized%20%3D%20R_x%20%2F%20R_x%5B0%5D%0A%0A%20%20%20%20print(f%22Autocorrelation%20at%20lag%200%3A%20%7BR_x%5B0%5D%3A.4f%7D%22)%0A%20%20%20%20print(f%22Autocorrelation%20at%20lag%201%3A%20%7BR_x%5B1%5D%3A.4f%7D%20(normalized%3A%20%7BR_x_normalized%5B1%5D%3A.4f%7D)%22)%0A%20%20%20%20for%20_lag%20in%20(%5B12%2C%2024%5D%20if%20is_monthly_data%20else%20%5B7%2C%2030%5D)%3A%0A%20%20%20%20%20%20%20%20print(f%22Autocorrelation%20at%20lag%20%7B_lag%7D%3A%20%7BR_x%5B_lag%5D%3A.4f%7D%20(normalized%3A%20%7BR_x_normalized%5B_lag%5D%3A.4f%7D)%22)%0A%0A%20%20%20%20%23%20Plot%20autocorrelation%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20plt.stem(lags%2C%20R_x_normalized%2C%20basefmt%3D'%20')%0A%20%20%20%20lag_label%20%3D%20'Lag%20%E2%84%93%20(months)'%20if%20is_monthly_data%20else%20'Lag%20%E2%84%93%20(days)'%0A%20%20%20%20plt.xlabel(lag_label)%0A%20%20%20%20plt.ylabel('Normalized%20Autocorrelation%20R_x%5B%E2%84%93%5D%20%2F%20R_x%5B0%5D')%0A%20%20%20%20plt.title('Autocorrelation%20Function%20of%20Temperature%20Anomalies')%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20xlim_max%20%3D%20min(max_lag%2C%2024)%20if%20is_monthly_data%20else%20min(max_lag%2C%20200)%0A%20%20%20%20plt.xlim(%5B0%2C%20xlim_max%5D)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20R_x%2C%20estimate_autocorrelation%0A%0A%0A%40app.cell%0Adef%20_(R_x%2C%20dates%2C%20fft%2C%20np%2C%20pd%2C%20plt)%3A%0A%20%20%20%20def%20compute_psd_from_acf(R%3A%20np.ndarray%2C%20fft_length%3A%20int%20%3D%20None)%20-%3E%20tuple%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Compute%20PSD%20from%20autocorrelation%20via%20FFT.%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20-----------%0A%20%20%20%20%20%20%20%20R%20%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20Autocorrelation%20R%5B0..L%5D%0A%20%20%20%20%20%20%20%20fft_length%20%3A%20int%2C%20optional%0A%20%20%20%20%20%20%20%20%20%20%20%20FFT%20length%20(default%3A%20next%20power%20of%202%20%3E%3D%202*len(R))%0A%20%20%20%20%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20--------%0A%20%20%20%20%20%20%20%20tuple%20(frequencies%2C%20PSD)%0A%20%20%20%20%20%20%20%20%20%20%20%20Frequencies%20(in%20cycles%20per%20sample%20period)%20and%20PSD%20values%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20L%20%3D%20len(R)%20-%201%0A%0A%20%20%20%20%20%20%20%20if%20fft_length%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Use%20next%20power%20of%202%20%3E%3D%202*L%0A%20%20%20%20%20%20%20%20%20%20%20%20fft_length%20%3D%20int(2**np.ceil(np.log2(2*L)))%0A%0A%20%20%20%20%20%20%20%20%23%20Construct%20symmetric%20sequence%20R%5B-L..L%5D%0A%20%20%20%20%20%20%20%20R_full%20%3D%20np.zeros(fft_length%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20R_full%5B%3AL%2B1%5D%20%3D%20R%20%20%23%20R%5B0..L%5D%0A%20%20%20%20%20%20%20%20R_full%5Bfft_length-L%3A%5D%20%3D%20R%5B1%3A%5D%5B%3A%3A-1%5D%20%20%23%20R%5B-L..-1%5D%20(symmetric)%0A%0A%20%20%20%20%20%20%20%20%23%20Compute%20PSD%20via%20FFT%0A%20%20%20%20%20%20%20%20S%20%3D%20np.real(fft.fft(R_full))%0A%0A%20%20%20%20%20%20%20%20%23%20Frequency%20axis%20(cycles%20per%20day)%0A%20%20%20%20%20%20%20%20freqs%20%3D%20fft.fftfreq(fft_length%2C%20d%3D1.0)%0A%0A%20%20%20%20%20%20%20%20return%20freqs%2C%20S%0A%0A%0A%20%20%20%20%23%20Compute%20PSD%0A%20%20%20%20K%20%3D%204096%20%20%23%20FFT%20length%20(power%20of%202)%0A%20%20%20%20freqs%2C%20S_x%20%3D%20compute_psd_from_acf(R_x%2C%20fft_length%3DK)%0A%0A%20%20%20%20%23%20Keep%20only%20nonnegative%20frequencies%0A%20%20%20%20nonneg_mask%20%3D%20freqs%20%3E%3D%200%0A%20%20%20%20freqs_pos%20%3D%20freqs%5Bnonneg_mask%5D%0A%20%20%20%20S_x_pos%20%3D%20S_x%5Bnonneg_mask%5D%0A%0A%20%20%20%20%23%20Determine%20frequency%20units%0A%20%20%20%20freq_check_psd%20%3D%20pd.infer_freq(dates)%0A%20%20%20%20median_days_psd%20%3D%20dates.to_series().diff().dt.days.median()%0A%20%20%20%20is_monthly_psd%20%3D%20bool((freq_check_psd%20and%20('M'%20in%20freq_check_psd%20or%20'MS'%20in%20freq_check_psd))%20or%20((dates.day%20%3D%3D%201).all()%20and%20median_days_psd%20%3E%3D%2028))%0A%20%20%20%20freq_unit%20%3D%20'cycles%2Fmonth'%20if%20is_monthly_psd%20else%20'cycles%2Fday'%0A%0A%20%20%20%20print(f%22PSD%20computed%20with%20FFT%20length%20%7BK%7D%22)%0A%20%20%20%20print(f%22Frequency%20resolution%3A%20%7Bfreqs%5B1%5D%3A.6f%7D%20%7Bfreq_unit%7D%22)%0A%20%20%20%20print(f%22Max%20frequency%3A%20%7Bfreqs_pos%5B-1%5D%3A.2f%7D%20%7Bfreq_unit%7D%20(Nyquist)%22)%0A%0A%20%20%20%20%23%20Plot%20PSD%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20plt.plot(freqs_pos%2C%20S_x_pos%2C%20'b-'%2C%20linewidth%3D1.5)%0A%20%20%20%20plt.xlabel(f'Frequency%20(%7Bfreq_unit%7D)')%0A%20%20%20%20plt.ylabel('Power%20Spectral%20Density')%0A%20%20%20%20plt.title('Power%20Spectral%20Density%20of%20Temperature%20Anomalies%20(from%20ACF)')%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20xlim_max_psd%20%3D%200.5%20if%20is_monthly_psd%20else%200.1%20%20%23%20Higher%20max%20for%20monthly%20(cycles%2Fmonth)%0A%20%20%20%20plt.xlim(%5B0%2C%20xlim_max_psd%5D)%20%20%23%20Focus%20on%20low%20frequencies%20(long-term%20patterns)%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20S_x%2C%0A%20%20%20%20%20%20%20%20S_x_pos%2C%0A%20%20%20%20%20%20%20%20compute_psd_from_acf%2C%0A%20%20%20%20%20%20%20%20freq_unit%2C%0A%20%20%20%20%20%20%20%20freqs_pos%2C%0A%20%20%20%20%20%20%20%20nonneg_mask%2C%0A%20%20%20%20%20%20%20%20xlim_max_psd%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20S_x%2C%0A%20%20%20%20S_x_pos%2C%0A%20%20%20%20fft%2C%0A%20%20%20%20freq_unit%2C%0A%20%20%20%20freqs_pos%2C%0A%20%20%20%20nonneg_mask%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20xlim_max_psd%2C%0A)%3A%0A%20%20%20%20def%20cepstral_minimum_phase_factor(S%3A%20np.ndarray%2C%20eps%3A%20float%20%3D%201e-10)%20-%3E%20np.ndarray%3A%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20Compute%20minimum-phase%20spectral%20factor%20via%20cepstral%20analysis.%0A%0A%20%20%20%20%20%20%20%20Parameters%3A%0A%20%20%20%20%20%20%20%20-----------%0A%20%20%20%20%20%20%20%20S%20%3A%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20Power%20spectral%20density%20(full%20symmetric%20spectrum%20from%20FFT)%0A%20%20%20%20%20%20%20%20eps%20%3A%20float%0A%20%20%20%20%20%20%20%20%20%20%20%20Small%20positive%20value%20to%20ensure%20S%20%3E%200%0A%20%20%20%20%0A%20%20%20%20%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20--------%0A%20%20%20%20%20%20%20%20np.ndarray%0A%20%20%20%20%20%20%20%20%20%20%20%20Complex%20frequency%20response%20H(%CF%89)%20(minimum-phase)%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20K%20%3D%20len(S)%0A%0A%20%20%20%20%20%20%20%20%23%20Ensure%20PSD%20is%20strictly%20positive%0A%20%20%20%20%20%20%20%20S_pos%20%3D%20np.maximum(S%2C%20eps)%0A%0A%20%20%20%20%20%20%20%20%23%20Log%20magnitude%20of%20spectral%20factor%3A%20log%20%7CH%7C%20%3D%20(1%2F2)%20log%20S%0A%20%20%20%20%20%20%20%20%23%20Since%20S%20is%20already%20symmetric%20(FFT%20of%20real%20autocorrelation)%2C%20log_mag%20is%20also%20symmetric%0A%20%20%20%20%20%20%20%20log_mag%20%3D%200.5%20*%20np.log(S_pos)%0A%0A%20%20%20%20%20%20%20%20%23%20Compute%20real%20cepstrum%3A%20c%20%3D%20IFFT(log%20%7CH%7C)%0A%20%20%20%20%20%20%20%20%23%20Since%20log_mag%20is%20symmetric%20(real)%2C%20the%20IFFT%20gives%20a%20real%2C%20symmetric%20cepstrum%0A%20%20%20%20%20%20%20%20c%20%3D%20np.real(fft.ifft(log_mag))%0A%0A%20%20%20%20%20%20%20%20%23%20Construct%20minimum-phase%20cepstrum%0A%20%20%20%20%20%20%20%20%23%20For%20minimum-phase%2C%20we%20keep%20only%20the%20causal%20part%20and%20double%20it%0A%20%20%20%20%20%20%20%20c_min%20%3D%20np.zeros(K%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20c_min%5B0%5D%20%3D%20c%5B0%5D%20%20%23%20DC%20component%0A%0A%20%20%20%20%20%20%20%20%23%20Causal%20part%3A%20double%20the%20positive-time%20cepstrum%0A%20%20%20%20%20%20%20%20%23%20Note%3A%20for%20a%20symmetric%20cepstrum%2C%20c%5Bk%5D%20%3D%20c%5B-k%5D%2C%20so%20we%20double%20the%20first%20half%0A%20%20%20%20%20%20%20%20if%20K%20%25%202%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Even%20length%3A%20indices%200%2C%201%2C%20...%2C%20K%2F2-1%2C%20K%2F2%0A%20%20%20%20%20%20%20%20%20%20%20%20c_min%5B1%3AK%2F%2F2%5D%20%3D%202%20*%20c%5B1%3AK%2F%2F2%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20c_min%5BK%2F%2F2%5D%20%3D%20c%5BK%2F%2F2%5D%20%20%23%20Nyquist%20frequency%20(not%20doubled)%0A%20%20%20%20%20%20%20%20%20%20%20%20c_min%5BK%2F%2F2%2B1%3A%5D%20%3D%200%20%20%23%20Anti-causal%20part%20set%20to%200%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Odd%20length%3A%20indices%200%2C%201%2C%20...%2C%20(K-1)%2F2%0A%20%20%20%20%20%20%20%20%20%20%20%20c_min%5B1%3A(K%2B1)%2F%2F2%5D%20%3D%202%20*%20c%5B1%3A(K%2B1)%2F%2F2%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20c_min%5B(K%2B1)%2F%2F2%3A%5D%20%3D%200%0A%0A%20%20%20%20%20%20%20%20%23%20Compute%20minimum-phase%20frequency%20response%3A%20H%20%3D%20exp(FFT(c_min))%0A%20%20%20%20%20%20%20%20log_H%20%3D%20fft.fft(c_min)%0A%20%20%20%20%20%20%20%20H%20%3D%20np.exp(log_H)%0A%0A%20%20%20%20%20%20%20%20return%20H%0A%0A%0A%20%20%20%20%23%20Apply%20cepstral%20factorization%0A%20%20%20%20%23%20S_x%20from%20compute_psd_from_acf%20is%20already%20the%20full%20symmetric%20spectrum%0A%20%20%20%20%23%20(FFT%20of%20symmetric%20autocorrelation%20gives%20symmetric%20PSD)%0A%20%20%20%20H%20%3D%20cepstral_minimum_phase_factor(S_x)%0A%0A%20%20%20%20%23%20Extract%20magnitude%20and%20use%20the%20same%20indexing%20as%20freqs_pos%0A%20%20%20%20H_mag%20%3D%20np.abs(H)%0A%20%20%20%20%23%20Use%20the%20same%20nonneg_mask%20that%20was%20used%20for%20freqs_pos%20to%20ensure%20alignment%0A%20%20%20%20H_mag_pos%20%3D%20H_mag%5Bnonneg_mask%5D%0A%0A%20%20%20%20%23%20Estimate%20noise%20variance%3A%20%CF%83%C2%B2%20%3D%20mean(S_x%20%2F%20%7CH%7C%C2%B2)%0A%20%20%20%20H_mag_sq%20%3D%20H_mag**2%0A%20%20%20%20sigma_sq%20%3D%20np.mean(S_x%20%2F%20(H_mag_sq%20%2B%201e-10))%0A%20%20%20%20sigma%20%3D%20np.sqrt(sigma_sq)%0A%0A%20%20%20%20print(f%22Estimated%20noise%20variance%20%CF%83%C2%B2%20%3D%20%7Bsigma_sq%3A.4f%7D%22)%0A%20%20%20%20print(f%22Estimated%20noise%20std%20%CF%83%20%3D%20%7Bsigma%3A.4f%7D%22)%0A%0A%20%20%20%20%23%20Plot%3A%20compare%20%7CH%7C%C2%B2%20to%20original%20PSD%0A%20%20%20%20%23%20Use%20the%20same%20frequency%20unit%20and%20limits%20as%20defined%20in%20the%20PSD%20cell%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20plt.plot(freqs_pos%2C%20S_x_pos%2C%20'b-'%2C%20linewidth%3D2%2C%20label%3D'Original%20PSD%20S_x(%CF%89)'%2C%20alpha%3D0.7)%0A%20%20%20%20plt.plot(freqs_pos%2C%20sigma_sq%20*%20H_mag_pos**2%2C%20'r--'%2C%20linewidth%3D2%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df'%CF%83%C2%B2%7CH(%CF%89)%7C%C2%B2%20(%CF%83%C2%B2%3D%7Bsigma_sq%3A.4f%7D)'%2C%20alpha%3D0.7)%0A%20%20%20%20plt.xlabel(f'Frequency%20(%7Bfreq_unit%7D)')%0A%20%20%20%20plt.ylabel('Power%20Spectral%20Density')%0A%20%20%20%20plt.title('Comparison%3A%20Original%20PSD%20vs.%20Minimum-Phase%20Factor%20%7CH(%CF%89)%7C%C2%B2')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlim(%5B0%2C%20xlim_max_psd%5D)%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%0A%20%20%20%20%23%20Plot%20magnitude%20of%20H%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20plt.plot(freqs_pos%2C%20H_mag_pos%2C%20'g-'%2C%20linewidth%3D2)%0A%20%20%20%20plt.xlabel(f'Frequency%20(%7Bfreq_unit%7D)')%0A%20%20%20%20plt.ylabel('%7CH(%CF%89)%7C')%0A%20%20%20%20plt.title('Magnitude%20of%20Minimum-Phase%20Spectral%20Factor%20H(%CF%89)')%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlim(%5B0%2C%20xlim_max_psd%5D)%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20(H%2C)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20H%2C%0A%20%20%20%20N%2C%0A%20%20%20%20anomalies%2C%0A%20%20%20%20compute_psd_from_acf%2C%0A%20%20%20%20dates%2C%0A%20%20%20%20estimate_autocorrelation%2C%0A%20%20%20%20fft%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20signal%2C%0A)%3A%0A%20%20%20%20%23%20Innovation%20filter%20in%20frequency%20domain%3A%20G%20%3D%201%2FH%0A%20%20%20%20G%20%3D%201.0%20%2F%20(H%20%2B%201e-10)%20%20%23%20Add%20small%20epsilon%20to%20avoid%20division%20by%20zero%0A%20%20%20%20g_full%20%3D%20np.real(fft.ifft(G))%0A%20%20%20%20%23%20Get%20time-domain%20impulse%20response%20via%20IFFT%0A%20%20%20%20filter_length%20%3D%20200%0A%20%20%20%20g%20%3D%20g_full%5B%3Afilter_length%5D%0A%20%20%20%20%23%20Truncate%20to%20manageable%20length%20(FIR%20filter)%0A%20%20%20%20g%20%3D%20g%20%2F%20np.sqrt(np.sum(g%20**%202))%0A%20%20%20%20print(f'Innovation%20filter%20length%3A%20%7Blen(g)%7D%20taps')%0A%20%20%20%20print(f'Filter%20energy%3A%20%7Bnp.sum(g%20**%202)%3A.4f%7D')%0A%20%20%20%20%23%20Normalize%20so%20that%20innovations%20have%20approximately%20unit%20variance%0A%20%20%20%20%23%20We%20can%20estimate%20this%20by%20applying%20to%20a%20test%20signal%20or%20by%20computing%20the%20gain%0A%20%20%20%20%23%20For%20now%2C%20normalize%20by%20the%20energy%20of%20the%20filter%0A%20%20%20%20innovations%20%3D%20signal.lfilter(g%2C%20%5B1.0%5D%2C%20anomalies)%0A%20%20%20%20innovations%20%3D%20innovations%20%2F%20np.std(innovations)%0A%20%20%20%20print(f'Innovations%20statistics%3A')%0A%20%20%20%20print(f'%20%20Mean%3A%20%7Bnp.mean(innovations)%3A.4f%7D')%0A%20%20%20%20print(f'%20%20Std%3A%20%7Bnp.std(innovations)%3A.4f%7D')%0A%20%20%20%20%23%20Apply%20innovation%20filter%20to%20anomalies%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(14%2C%208)%2C%20sharex%3DTrue)%0A%20%20%20%20window_days%20%3D%20180%0A%20%20%20%20%23%20Normalize%20innovations%20to%20have%20unit%20variance%20(for%20comparison)%0A%20%20%20%20_start_idx%20%3D%20max(0%2C%20N%20-%20window_days)%0A%20%20%20%20_axes%5B0%5D.plot(dates%5B_start_idx%3A%5D%2C%20anomalies%5B_start_idx%3A%5D%2C%20'b-'%2C%20alpha%3D0.7%2C%20linewidth%3D0.8)%0A%20%20%20%20_axes%5B0%5D.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20alpha%3D0.3)%0A%20%20%20%20_axes%5B0%5D.set_ylabel('Temperature%20Anomaly%20(%C2%B0C)')%0A%20%20%20%20_axes%5B0%5D.set_title('Temperature%20Anomalies')%0A%20%20%20%20_axes%5B0%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20%23%20Plot%3A%20anomalies%20vs%20innovations%20(short%20segment)%0A%20%20%20%20_axes%5B1%5D.plot(dates%5B_start_idx%3A%5D%2C%20innovations%5B_start_idx%3A%5D%2C%20'r-'%2C%20alpha%3D0.7%2C%20linewidth%3D0.8)%0A%20%20%20%20_axes%5B1%5D.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20alpha%3D0.3)%0A%20%20%20%20%23%20Show%206%20months%0A%20%20%20%20_axes%5B1%5D.set_ylabel('Innovations%20(normalized)')%0A%20%20%20%20_axes%5B1%5D.set_xlabel('Date')%0A%20%20%20%20_axes%5B1%5D.set_title('Innovations%20(Whitened%20Signal)')%0A%20%20%20%20_axes%5B1%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20max_lag_innov%20%3D%2050%0A%20%20%20%20R_innov%20%3D%20estimate_autocorrelation(innovations%2C%20max_lag_innov)%0A%20%20%20%20R_innov_normalized%20%3D%20R_innov%20%2F%20R_innov%5B0%5D%20if%20R_innov%5B0%5D%20%3E%200%20else%20R_innov%0A%20%20%20%20print(f'Innovation%20ACF%20at%20lag%200%3A%20%7BR_innov%5B0%5D%3A.4f%7D')%0A%20%20%20%20print(f'Innovation%20ACF%20at%20lag%201%3A%20%7BR_innov%5B1%5D%3A.4f%7D%20(normalized%3A%20%7BR_innov_normalized%5B1%5D%3A.4f%7D)')%0A%20%20%20%20print(f'Innovation%20ACF%20at%20lag%205%3A%20%7BR_innov%5B5%5D%3A.4f%7D%20(normalized%3A%20%7BR_innov_normalized%5B5%5D%3A.4f%7D)')%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20lags_innov%20%3D%20np.arange(len(R_innov))%0A%20%20%20%20plt.stem(lags_innov%2C%20R_innov_normalized%2C%20basefmt%3D'%20')%0A%20%20%20%20plt.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'-'%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlabel('Lag%20%E2%84%93%20(days)')%0A%20%20%20%20plt.ylabel('Normalized%20Autocorrelation')%0A%20%20%20%20plt.title('Autocorrelation%20Function%20of%20Innovations%20(should%20be%20%E2%89%88%20%CE%B4%5B%E2%84%93%5D)')%0A%20%20%20%20%23%20Check%20whiteness%3A%20estimate%20autocorrelation%20of%20innovations%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlim(%5B0%2C%20max_lag_innov%5D)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20freqs_innov%2C%20S_innov%20%3D%20compute_psd_from_acf(R_innov%2C%20fft_length%3D1024)%0A%20%20%20%20nonneg_mask_innov%20%3D%20freqs_innov%20%3E%3D%200%0A%20%20%20%20freqs_innov_pos%20%3D%20freqs_innov%5Bnonneg_mask_innov%5D%0A%20%20%20%20S_innov_pos%20%3D%20S_innov%5Bnonneg_mask_innov%5D%0A%20%20%20%20%23%20Plot%20ACF%20of%20innovations%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20plt.plot(freqs_innov_pos%2C%20S_innov_pos%2C%20'r-'%2C%20linewidth%3D2)%0A%20%20%20%20plt.axhline(y%3Dnp.mean(S_innov_pos)%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20label%3Df'Mean%20PSD%20%3D%20%7Bnp.mean(S_innov_pos)%3A.4f%7D'%2C%20alpha%3D0.7)%0A%20%20%20%20plt.xlabel('Frequency%20(cycles%2Fday)')%0A%20%20%20%20plt.ylabel('Power%20Spectral%20Density')%0A%20%20%20%20plt.title('PSD%20of%20Innovations%20(should%20be%20relatively%20flat%20for%20white%20noise)')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlim(%5B0%2C%200.5%5D)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20%23%20Estimate%20PSD%20of%20innovations%0A%20%20%20%20%23%20Plot%20PSD%20of%20innovations%0A%20%20%20%20plt.show()%0A%20%20%20%20return%20R_innov_normalized%2C%20innovations%2C%20lags_innov%2C%20max_lag_innov%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20AutoReg%2C%0A%20%20%20%20HAS_STATSMODELS%2C%0A%20%20%20%20N%2C%0A%20%20%20%20R_innov_normalized%2C%0A%20%20%20%20R_x%2C%0A%20%20%20%20anomalies%2C%0A%20%20%20%20dates%2C%0A%20%20%20%20estimate_autocorrelation%2C%0A%20%20%20%20innovations%2C%0A%20%20%20%20lags_innov%2C%0A%20%20%20%20max_lag_innov%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A)%3A%0A%20%20%20%20%23%20Choose%20AR%20order%0A%20%20%20%20ar_order%20%3D%2010%0A%20%20%20%20if%20HAS_STATSMODELS%3A%0A%20%20%20%20%20%20%20%20print(f'Fitting%20AR(%7Bar_order%7D)%20model%20using%20statsmodels...')%0A%20%20%20%20%20%20%20%20ar_model%20%3D%20AutoReg(anomalies%2C%20lags%3Dar_order%2C%20trend%3D'n')%20%20%23%20Option%20A%3A%20Use%20statsmodels%0A%20%20%20%20%20%20%20%20ar_result%20%3D%20ar_model.fit()%0A%20%20%20%20%20%20%20%20ar_coeffs%20%3D%20-ar_result.params.values%0A%20%20%20%20%20%20%20%20print(f'AR%20coefficients%3A%20%7Bar_coeffs%7D')%0A%20%20%20%20else%3A%20%20%23%20statsmodels%20uses%20opposite%20sign%20convention%0A%20%20%20%20%20%20%20%20print(f'Fitting%20AR(%7Bar_order%7D)%20model%20using%20Yule-Walker%20equations...')%0A%20%20%20%20%20%20%20%20R_matrix%20%3D%20np.zeros((ar_order%2C%20ar_order))%0A%20%20%20%20%20%20%20%20r_vector%20%3D%20np.zeros(ar_order)%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(ar_order)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20j%20in%20range(ar_order)%3A%20%20%23%20Option%20B%3A%20Solve%20Yule-Walker%20equations%20directly%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20lag%20%3D%20abs(i%20-%20j)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20lag%20%3C%20len(R_x)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20R_matrix%5Bi%2C%20j%5D%20%3D%20R_x%5Blag%5D%20%20%23%20Build%20Yule-Walker%20system%3A%20R%20a%20%3D%20-r%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20i%20%2B%201%20%3C%20len(R_x)%3A%20%20%23%20where%20R%20is%20Toeplitz%20matrix%20of%20autocorrelations%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20r_vector%5Bi%5D%20%3D%20R_x%5Bi%20%2B%201%5D%0A%20%20%20%20%20%20%20%20ar_coeffs%20%3D%20-np.linalg.solve(R_matrix%2C%20r_vector)%0A%20%20%20%20%20%20%20%20print(f'AR%20coefficients%3A%20%7Bar_coeffs%7D')%0A%20%20%20%20predictions%20%3D%20np.zeros_like(anomalies)%0A%20%20%20%20predictions%5B%3Aar_order%5D%20%3D%20anomalies%5B%3Aar_order%5D%0A%20%20%20%20for%20n%20in%20range(ar_order%2C%20len(anomalies))%3A%0A%20%20%20%20%20%20%20%20pred%20%3D%200.0%0A%20%20%20%20%20%20%20%20for%20k%20in%20range(ar_order)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pred%20-%3D%20ar_coeffs%5Bk%5D%20*%20anomalies%5Bn%20-%20k%20-%201%5D%0A%20%20%20%20%20%20%20%20predictions%5Bn%5D%20%3D%20pred%0A%20%20%20%20pred_errors%20%3D%20anomalies%20-%20predictions%0A%20%20%20%20pred_errors_normalized%20%3D%20pred_errors%20%2F%20np.std(pred_errors)%20%20%23%20Solve%3A%20R%20a%20%3D%20-r%20%20%3D%3E%20%20a%20%3D%20-R%5E%7B-1%7D%20r%0A%20%20%20%20print(f'%5CnPrediction%20error%20statistics%3A')%0A%20%20%20%20print(f'%20%20Mean%3A%20%7Bnp.mean(pred_errors)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20print(f'%20%20Std%3A%20%7Bnp.std(pred_errors)%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20print(f'%20%20RMSE%3A%20%7Bnp.sqrt(np.mean(pred_errors%20**%202))%3A.4f%7D%20%C2%B0C')%0A%20%20%20%20%23%20Build%20one-step-ahead%20predictions%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(14%2C%208)%2C%20sharex%3DTrue)%0A%20%20%20%20window_days_pred%20%3D%2090%20%20%23%20Initialize%20with%20actual%20values%0A%20%20%20%20start_idx_pred%20%3D%20max(ar_order%2C%20N%20-%20window_days_pred)%0A%20%20%20%20_axes%5B0%5D.plot(dates%5Bstart_idx_pred%3A%5D%2C%20anomalies%5Bstart_idx_pred%3A%5D%2C%20'b-'%2C%20alpha%3D0.7%2C%20linewidth%3D1.5%2C%20label%3D'Anomalies')%0A%20%20%20%20_axes%5B0%5D.plot(dates%5Bstart_idx_pred%3A%5D%2C%20predictions%5Bstart_idx_pred%3A%5D%2C%20'r--'%2C%20alpha%3D0.8%2C%20linewidth%3D1.5%2C%20label%3D'One-step-ahead%20prediction')%20%20%23%20x%CC%82%5Bn%7Cn-1%5D%20%3D%20-%CE%A3%20a_k%20x%5Bn-k%5D%0A%20%20%20%20_axes%5B0%5D.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20alpha%3D0.3)%0A%20%20%20%20_axes%5B0%5D.set_ylabel('Temperature%20Anomaly%20(%C2%B0C)')%0A%20%20%20%20_axes%5B0%5D.set_title(f'AR(%7Bar_order%7D)%20One-Step-Ahead%20Prediction')%0A%20%20%20%20_axes%5B0%5D.legend()%0A%20%20%20%20_axes%5B0%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20%23%20Prediction%20errors%0A%20%20%20%20_axes%5B1%5D.plot(dates%5Bstart_idx_pred%3A%5D%2C%20pred_errors%5Bstart_idx_pred%3A%5D%2C%20'g-'%2C%20alpha%3D0.7%2C%20linewidth%3D0.8)%0A%20%20%20%20_axes%5B1%5D.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'--'%2C%20alpha%3D0.3)%0A%20%20%20%20%23%20Normalize%20prediction%20errors%20for%20comparison%0A%20%20%20%20_axes%5B1%5D.set_ylabel('Prediction%20Error%20(%C2%B0C)')%0A%20%20%20%20_axes%5B1%5D.set_xlabel('Date')%0A%20%20%20%20_axes%5B1%5D.set_title('Prediction%20Errors')%0A%20%20%20%20_axes%5B1%5D.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20%23%20Plot%3A%20anomalies%20vs%20predictions%20(short%20window)%0A%20%20%20%20plt.hist(innovations%2C%20bins%3D50%2C%20alpha%3D0.5%2C%20label%3D'Innovations%20(cepstral)'%2C%20density%3DTrue)%0A%20%20%20%20plt.hist(pred_errors_normalized%2C%20bins%3D50%2C%20alpha%3D0.5%2C%20label%3D'Prediction%20errors%20(AR)'%2C%20density%3DTrue)%0A%20%20%20%20%23%20Show%203%20months%0A%20%20%20%20plt.xlabel('Normalized%20value')%0A%20%20%20%20plt.ylabel('Density')%0A%20%20%20%20plt.title('Comparison%3A%20Innovations%20vs.%20AR%20Prediction%20Errors')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20plt.show()%0A%20%20%20%20R_pred%20%3D%20estimate_autocorrelation(pred_errors%2C%20max_lag_innov)%0A%20%20%20%20R_pred_normalized%20%3D%20R_pred%20%2F%20R_pred%5B0%5D%20if%20R_pred%5B0%5D%20%3E%200%20else%20R_pred%0A%20%20%20%20plt.figure(figsize%3D(12%2C%206))%0A%20%20%20%20lags_pred%20%3D%20np.arange(len(R_pred))%0A%20%20%20%20plt.stem(lags_pred%2C%20R_pred_normalized%2C%20basefmt%3D'%20'%2C%20label%3D'Prediction%20errors')%0A%20%20%20%20plt.stem(lags_innov%2C%20R_innov_normalized%2C%20basefmt%3D'%20'%2C%20markerfmt%3D'ro'%2C%20linefmt%3D'r--'%2C%20label%3D'Innovations')%0A%20%20%20%20plt.axhline(y%3D0%2C%20color%3D'k'%2C%20linestyle%3D'-'%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlabel('Lag%20%E2%84%93%20(days)')%0A%20%20%20%20plt.ylabel('Normalized%20Autocorrelation')%0A%20%20%20%20plt.title('ACF%20Comparison%3A%20Prediction%20Errors%20vs.%20Innovations')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20plt.xlim(%5B0%2C%20max_lag_innov%5D)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20%23%20Compare%20prediction%20errors%20to%20innovations%0A%20%20%20%20%23%20ACF%20of%20prediction%20errors%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Interpretation%20and%20Discussion%0A%0A%20%20%20%20%23%23%23%20Signal%20Processing%20Perspective%0A%0A%20%20%20%20**Temperature%20anomalies%20as%20filtered%20white%20noise%3A**%0A%20%20%20%20The%20temperature%20anomalies%20x%5Bn%5D%20can%20be%20viewed%20as%20the%20output%20of%20a%20linear%20time-invariant%20(LTI)%0A%20%20%20%20system%20H(z)%20driven%20by%20white%20noise%20w%5Bn%5D%20(weather%20shocks)%3A%0A%0A%20%20%20%20%20%20%20%20x%5Bn%5D%20%3D%20H(z)%20*%20w%5Bn%5D%20%20%20%20%E2%9F%BA%20%20%20%20X(%CF%89)%20%3D%20H(%CF%89)%20W(%CF%89)%0A%0A%20%20%20%20The%20power%20spectral%20density%20relationship%20is%3A%0A%0A%20%20%20%20%20%20%20%20S_x(%CF%89)%20%3D%20%CF%83%C2%B2%20%7CH(%CF%89)%7C%C2%B2%0A%0A%20%20%20%20where%20%CF%83%C2%B2%20is%20the%20variance%20of%20the%20white%20noise%20input.%0A%0A%20%20%20%20**Cepstral%20spectral%20factorization%3A**%0A%20%20%20%20The%20cepstrum%20(inverse%20FFT%20of%20log-magnitude%20spectrum)%20allows%20us%20to%20decompose%20the%20log-spectrum%0A%20%20%20%20into%20causal%20and%20anti-causal%20components.%20By%20doubling%20the%20causal%20part%20and%20zeroing%20the%20anti-causal%0A%20%20%20%20part%2C%20we%20obtain%20a%20minimum-phase%20spectral%20factor%20H(%CF%89)%20such%20that%20%7CH(%CF%89)%7C%C2%B2%20matches%20S_x(%CF%89).%0A%0A%20%20%20%20**Innovation%20(whitening)%20filter%3A**%0A%20%20%20%20The%20innovation%20filter%20G(%CF%89)%20%3D%201%2FH(%CF%89)%20inverts%20the%20shaping%20filter%2C%20extracting%20the%20white%20innovations%3A%0A%0A%20%20%20%20%20%20%20%20w%5Bn%5D%20%3D%20G(z)%20*%20x%5Bn%5D%20%20%20%20%E2%9F%BA%20%20%20%20W(%CF%89)%20%3D%20G(%CF%89)%20X(%CF%89)%20%3D%20X(%CF%89)%20%2F%20H(%CF%89)%0A%0A%20%20%20%20The%20innovations%20represent%20the%20unpredictable%20component%E2%80%94the%20%22surprise%22%20at%20each%20time%20step%20that%0A%20%20%20%20cannot%20be%20predicted%20from%20past%20values.%0A%0A%20%20%20%20**AR%20prediction%20and%20the%20innovation%20representation%3A**%0A%20%20%20%20An%20AR(p)%20model%20represents%20the%20anomalies%20as%3A%0A%0A%20%20%20%20%20%20%20%20x%5Bn%5D%20%3D%20-a%E2%82%81x%5Bn-1%5D%20-%20a%E2%82%82x%5Bn-2%5D%20-%20...%20-%20a%E2%82%9Ax%5Bn-p%5D%20%2B%20w%5Bn%5D%0A%0A%20%20%20%20The%20one-step-ahead%20predictor%20is%3A%0A%0A%20%20%20%20%20%20%20%20x%CC%82%5Bn%7Cn-1%5D%20%3D%20-a%E2%82%81x%5Bn-1%5D%20-%20a%E2%82%82x%5Bn-2%5D%20-%20...%20-%20a%E2%82%9Ax%5Bn-p%5D%0A%0A%20%20%20%20The%20prediction%20error%20e%5Bn%5D%20%3D%20x%5Bn%5D%20-%20x%CC%82%5Bn%7Cn-1%5D%20equals%20the%20innovation%20w%5Bn%5D.%20For%20the%20optimal%0A%20%20%20%20linear%20predictor%2C%20the%20prediction%20error%20is%20white%20(uncorrelated)%2C%20which%20is%20exactly%20what%20we%0A%20%20%20%20observe%3A%20the%20innovations%20from%20the%20cepstral%20method%20and%20the%20AR%20prediction%20errors%20are%20both%0A%20%20%20%20approximately%20white.%0A%0A%20%20%20%20**Key%20relationships%3A**%0A%0A%20%20%20%201.%20**Autocorrelation%20%E2%86%94%20PSD%3A**%20R_x%5B%E2%84%93%5D%20%E2%9F%B7%20S_x(%CF%89)%20via%20Fourier%20transform%0A%20%20%20%202.%20**PSD%20%E2%86%94%20Spectral%20factor%3A**%20S_x(%CF%89)%20%3D%20%CF%83%C2%B2%7CH(%CF%89)%7C%C2%B2%20via%20cepstral%20factorization%0A%20%20%20%203.%20**Spectral%20factor%20%E2%86%94%20Innovation%20filter%3A**%20G(%CF%89)%20%3D%201%2FH(%CF%89)%20(inverse%20system)%0A%20%20%20%204.%20**AR%20model%20%E2%86%94%20Innovation%3A**%20The%20AR%20prediction%20error%20equals%20the%20innovation%0A%0A%20%20%20%20**Physical%20interpretation%3A**%0A%0A%20%20%20%20-%20The%20autocorrelation%20structure%20captures%20persistent%20weather%20patterns%20(e.g.%2C%20multi-day%0A%20%20%20%20%20%20weather%20systems%2C%20seasonal%20persistence)%0A%20%20%20%20-%20The%20innovations%20represent%20day-to-day%20weather%20%22shocks%22%20that%20are%20unpredictable%20from%0A%20%20%20%20%20%20past%20observations%0A%20%20%20%20-%20The%20shaping%20filter%20H(%CF%89)%20describes%20how%20these%20shocks%20are%20%22colored%22%20by%20the%20climate%20system%0A%20%20%20%20%20%20to%20produce%20the%20observed%20temperature%20anomalies%0A%0A%20%20%20%20This%20framework%20connects%20classical%20signal%20processing%20(spectral%20factorization%2C%20whitening)%0A%20%20%20%20with%20time-series%20analysis%20(AR%20modeling%2C%20prediction)%20and%20provides%20a%20unified%20view%20of%0A%20%20%20%20temperature%20variability%20as%20filtered%20white%20noise.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20print(%22%5Cn%22%20%2B%20%22%3D%22*70)%0A%20%20%20%20print(%22Analysis%20complete!%20See%20plots%20above%20for%20visualization%20of%20each%20step.%22)%0A%20%20%20%20print(%22%3D%22*70)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A