# Fourier Regression and Seasonal Decomposition of Global Mean Sea Level
Related notes:
- [[MATH310S26-Day7-Notes]] - Introduction to Fourier series
- [[MATH310S26-Day10-Notes]] - Complex Fourier series
- [[MATH310S26-Day11-Notes]] - Fourier regression theory
- [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)]] - Normal equations template
- [[MATH310S26-Day6-Work]] - PCA code template
## Day 11 recap (see [[MATH310S26-Day11-Notes]])
- Formulated [Fourier series](https://en.wikipedia.org/wiki/Fourier_series) as a regression problem with [trigonometric basis functions](https://en.wikipedia.org/wiki/Trigonometric_functions)
- Applied to [Global Mean Sea Level (GMSL)](https://climate.nasa.gov/vital-signs/sea-level/) data to separate seasonality from trend
- Found that trend parameters (~3.3 mm/year) are essentially identical with or without seasonal removal
- Discovered [El Niño](https://en.wikipedia.org/wiki/El_Niño–Southern_Oscillation) patterns in residuals after removing annual/semi-annual cycles
- Showed that [R²](https://en.wikipedia.org/wiki/Coefficient_of_determination) improves significantly when fitting trends to de-seasonalized data
## Day 13 goals
- Implement Fourier regression as [ordinary least squares](https://en.wikipedia.org/wiki/Ordinary_least_squares) with [design matrices](https://en.wikipedia.org/wiki/Design_matrix)
- Extract seasonal components (annual and semi-annual) from GMSL [time series](https://en.wikipedia.org/wiki/Time_series)
- Compare trend estimation with and without [seasonal adjustment](https://en.wikipedia.org/wiki/Seasonal_adjustment)
- Visualize decomposition and identify climate patterns in residuals
# From theory to practice: Fourier regression as OLS
When we have [time series](https://en.wikipedia.org/wiki/Time_series) data $(t_i, y_i)$ for $i = 1, \ldots, m$, we can model it as:
$y_i \approx a_0 + \sum_{k=1}^{N} \left[ a_k \cos(2\pi k f_0 t_i) + b_k \sin(2\pi k f_0 t_i) \right]$
This is **[linear regression](https://en.wikipedia.org/wiki/Linear_regression)** with trigonometric features! The [design matrix](https://en.wikipedia.org/wiki/Design_matrix) has columns:
- Column 1: ones (for $a_0$)
- Column 2: $\cos(2\pi f_0 t)$ values
- Column 3: $\sin(2\pi f_0 t)$ values
- Column 4: $\cos(4\pi f_0 t)$ values
- Column 5: $\sin(4\pi f_0 t)$ values
- ... and so on
For sea level data with $t$ in years, we choose:
- $f_0 = 1$ year⁻¹ (fundamental frequency)
- $N = 2$ (capture annual and semi-annual cycles)
# Global Mean Sea Level Analysis
We'll work with GMSL data from 1993-2024 (source: [NASA Sea Level Portal](https://sealevel.nasa.gov/), see: [global_mean_sea_level_1993-2024.csv](https://www.dropbox.com/scl/fi/9uqs42k9k3tgqwl45gdn9/global_mean_sea_level_1993-2024.csv?rlkey=d62q93mjrso72ql76dhzpolg4&dl=0), or Canvas for the data) and apply the [time series decomposition](https://en.wikipedia.org/wiki/Decomposition_of_time_series):
$\text{Data} = \text{Trend} + \text{Seasonality} + \text{Error}$
## Python (NumPy + Matplotlib)
Step 1. Load and explore the data
```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Load the CSV file
df = pd.read_csv("global_mean_sea_level_1993-2024.csv")
df = df.dropna(subset=['GMSLWithGIA']) # Remove any NaN values
print(f"Data shape: {df.shape}") # expect ~1169 rows
```
Purpose: Load GMSL data, removing any missing values.
Step 2. Extract time and sea level variables
```python
t = df['YearPlusFraction'].values
S = df['GMSLWithGIA'].values
t_centered = t - t.mean() # Center for numerical stability
print(f"Time range: {t.min():.1f} to {t.max():.1f} years")
```
Purpose: Prepare time (in years) and sea level (in mm) arrays. Centering improves numerical conditioning.
Step 3. Quick visualization of raw data
```python
plt.figure(figsize=(10, 4))
plt.plot(t, S, 'b-', alpha=0.7, linewidth=0.5)
plt.xlabel('Year'); plt.ylabel('Sea Level (mm)')
plt.title('Global Mean Sea Level (1993-2024)')
plt.grid(True, alpha=0.3); plt.show()
```
Purpose: Visualize the time series to see trend and oscillations.
Step 4. Build Fourier design matrix
```python
# Annual (k=1) and semi-annual (k=2) frequencies
X_fourier = np.ones((len(t), 5))
X_fourier[:, 1] = np.cos(2*np.pi*1*t) # cos(2πt)
X_fourier[:, 2] = np.sin(2*np.pi*1*t) # sin(2πt)
X_fourier[:, 3] = np.cos(2*np.pi*2*t) # cos(4πt)
X_fourier[:, 4] = np.sin(2*np.pi*2*t) # sin(4πt)
print(f"Design matrix shape: {X_fourier.shape}")
```
Purpose: Create design matrix with intercept and two harmonic pairs.
Step 5. Solve for Fourier coefficients
```python
# Normal equations: X'X β = X'S
beta_fourier = np.linalg.solve(X_fourier.T @ X_fourier, X_fourier.T @ S)
S_seasonal = X_fourier @ beta_fourier
print(f"Fourier coefficients: {beta_fourier.round(2)}")
```
Purpose: Find coefficients via [least squares](https://en.wikipedia.org/wiki/Least_squares). The seasonal component is the fitted values.
Step 6. Extract and visualize seasonal component
```python
plt.figure(figsize=(10, 4))
plt.plot(t[:100], S[:100], 'b-', alpha=0.7, label='Raw data')
plt.plot(t[:100], S_seasonal[:100], 'r-', linewidth=2, label='Seasonal fit')
plt.xlabel('Year'); plt.ylabel('Sea Level (mm)')
plt.title('First 100 points: Raw vs Seasonal Fit')
plt.legend(); plt.grid(True, alpha=0.3); plt.show()
```
Purpose: Zoom in to see how well the Fourier model captures oscillations.
Step 7. Compute residuals (de-seasonalized data)
```python
residuals = S - S_seasonal
print(f"Residual std before: {S.std():.2f} mm")
print(f"Residual std after: {residuals.std():.2f} mm")
```
Purpose: Remove seasonal component. Standard deviation should decrease.
Step 8. Fit linear trend to residuals
```python
# Design matrix for linear trend (using centered time)
X_linear = np.column_stack([np.ones_like(t_centered), t_centered])
beta_linear = np.linalg.solve(X_linear.T @ X_linear, X_linear.T @ residuals)
trend_linear = X_linear @ beta_linear
print(f"Sea level rise rate: {beta_linear[1]:.3f} mm/year")
```
Purpose: Estimate trend from de-seasonalized data. Expect ~3.3 mm/year.
Step 9. Fit quadratic trend to residuals
```python
# Design matrix for quadratic trend
X_quad = np.column_stack([np.ones_like(t_centered), t_centered, t_centered**2])
beta_quad = np.linalg.solve(X_quad.T @ X_quad, X_quad.T @ residuals)
trend_quad = X_quad @ beta_quad
print(f"Acceleration: {2*beta_quad[2]:.4f} mm/year²")
```
Purpose: Test for acceleration in sea level rise.
Step 10. Calculate R² for both models
```python
SS_tot = np.sum((residuals - residuals.mean())**2)
SS_res_linear = np.sum((residuals - trend_linear)**2)
SS_res_quad = np.sum((residuals - trend_quad)**2)
R2_linear = 1 - SS_res_linear/SS_tot
R2_quad = 1 - SS_res_quad/SS_tot
print(f"R² linear: {R2_linear:.4f}, R² quadratic: {R2_quad:.4f}")
```
Purpose: Quantify how well trends explain de-seasonalized data. Higher [R²](https://en.wikipedia.org/wiki/Coefficient_of_determination) indicates better fit.
Step 11. Compare with direct fit (no seasonal removal)
```python
# Fit directly to original data
beta_direct = np.linalg.solve(X_linear.T @ X_linear, X_linear.T @ S)
print(f"Direct rate: {beta_direct[1]:.3f} mm/year")
print(f"De-seasoned rate: {beta_linear[1]:.3f} mm/year")
```
Purpose: Show that trend estimates are nearly identical with or without seasonal removal.
Step 12. Visualize complete decomposition
```python
fig, axes = plt.subplots(3, 1, figsize=(12, 8))
# Original data
axes[0].plot(t, S, 'b-', alpha=0.7, linewidth=0.5)
axes[0].set_ylabel('Sea Level (mm)')
axes[0].set_title('Original Data')
axes[0].grid(True, alpha=0.3)
# Seasonal component
axes[1].plot(t, S_seasonal - S_seasonal.mean(), 'g-', linewidth=1)
axes[1].set_ylabel('Seasonal (mm)')
axes[1].set_title('Seasonal Component (mean-centered)')
axes[1].grid(True, alpha=0.3)
# Residuals with trend
axes[2].plot(t, residuals, 'gray', alpha=0.5, linewidth=0.5)
axes[2].plot(t, trend_quad, 'r-', linewidth=2, label='Quadratic trend')
axes[2].set_xlabel('Year'); axes[2].set_ylabel('Residual (mm)')
axes[2].set_title('De-seasonalized Data with Trend')
axes[2].legend(); axes[2].grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
```
Purpose: Show full decomposition: data = seasonal + trend + noise.
Step 13. Examine final residuals for patterns
```python
final_residuals = residuals - trend_quad
plt.figure(figsize=(10, 4))
plt.plot(t, final_residuals, 'k-', alpha=0.7, linewidth=0.5)
plt.axhline(0, color='r', linestyle='--', alpha=0.5)
plt.xlabel('Year'); plt.ylabel('Final Residual (mm)')
plt.title('Residuals after removing seasonal + quadratic trend')
plt.grid(True, alpha=0.3); plt.show()
# Look for patterns - El Niño events appear as peaks around 1997-98, 2015-16
# Note: While residuals show ENSO structure, the physics is more complex than
# simple seasonal removal. ENSO is quasi-periodic (2-7 years), not strictly annual,
# and interacts nonlinearly with other climate modes like the [Pacific Decadal Oscillation](https://en.wikipedia.org/wiki/Pacific_decadal_oscillation)
# and [Indian Ocean Dipole](https://en.wikipedia.org/wiki/Indian_Ocean_Dipole).
```
Purpose: Remaining structure reveals climate oscillations like El Niño, though capturing ENSO fully requires more sophisticated methods than Fourier decomposition.
## R Implementation
Step 1. Load and prepare data
```r
df <- read.csv("global_mean_sea_level_1993-2024.csv")
df <- df[!is.na(df$GMSLWithGIA), ]
t <- df$YearPlusFraction
S <- df$GMSLWithGIA
t_centered <- t - mean(t)
cat("Data points:", length(S), "\n")
```
Purpose: Load GMSL data and prepare variables.
Step 2. Build Fourier design matrix
```r
n <- length(t)
X_fourier <- matrix(1, nrow=n, ncol=5)
X_fourier[, 2] <- cos(2*pi*1*t) # Annual cosine
X_fourier[, 3] <- sin(2*pi*1*t) # Annual sine
X_fourier[, 4] <- cos(2*pi*2*t) # Semi-annual cosine
X_fourier[, 5] <- sin(2*pi*2*t) # Semi-annual sine
```
Purpose: Create design matrix with trigonometric basis functions.
Step 3. Solve for Fourier coefficients
```r
beta_fourier <- solve(t(X_fourier) %*% X_fourier) %*% (t(X_fourier) %*% S)
S_seasonal <- X_fourier %*% beta_fourier
cat("Fourier coefficients:\n")
print(round(beta_fourier, 2))
```
Purpose: Find seasonal component via [normal equations](https://en.wikipedia.org/wiki/Normal_equation).
Step 4. De-seasonalize and fit trends
```r
residuals <- S - S_seasonal
# Linear trend on centered time
X_linear <- cbind(1, t_centered)
beta_linear <- solve(t(X_linear) %*% X_linear) %*% (t(X_linear) %*% residuals)
trend_linear <- X_linear %*% beta_linear
cat("Sea level rise rate:", round(beta_linear[2], 3), "mm/year\n")
```
Purpose: Remove seasonality and estimate linear trend.
Step 5. Fit quadratic trend
```r
X_quad <- cbind(1, t_centered, t_centered^2)
beta_quad <- solve(t(X_quad) %*% X_quad) %*% (t(X_quad) %*% residuals)
trend_quad <- X_quad %*% beta_quad
cat("Acceleration:", round(2*beta_quad[3], 4), "mm/year²\n")
```
Purpose: Test for acceleration in sea level rise.
Step 6. Calculate R² values
```r
SS_tot <- sum((residuals - mean(residuals))^2)
SS_res_linear <- sum((residuals - trend_linear)^2)
SS_res_quad <- sum((residuals - trend_quad)^2)
R2_linear <- 1 - SS_res_linear/SS_tot
R2_quad <- 1 - SS_res_quad/SS_tot
cat("R² linear:", round(R2_linear, 4), "\n")
cat("R² quadratic:", round(R2_quad, 4), "\n")
```
Purpose: Quantify goodness of fit for trend models.
Step 7. Visualization
```r
par(mfrow=c(3,1), mar=c(4,4,2,1))
# Original data
plot(t, S, type='l', col='blue', main='Original GMSL Data',
xlab='', ylab='Sea Level (mm)', lwd=0.5)
grid()
# Seasonal component
plot(t, S_seasonal - mean(S_seasonal), type='l', col='green',
main='Seasonal Component', xlab='', ylab='Seasonal (mm)')
grid()
# De-seasonalized with trend
plot(t, residuals, type='l', col='gray', main='De-seasonalized with Trend',
xlab='Year', ylab='Residual (mm)', lwd=0.5)
lines(t, trend_quad, col='red', lwd=2)
grid()
```
Purpose: Visualize the complete decomposition.
Step 8. Final residual analysis
```r
final_residuals <- residuals - trend_quad
par(mfrow=c(1,1))
plot(t, final_residuals, type='l', main='Final Residuals',
xlab='Year', ylab='Residual (mm)', lwd=0.5)
abline(h=0, col='red', lty=2)
grid()
# Peaks around 1998, 2016 correspond to strong El Niño events
```
Purpose: Identify climate patterns in residuals.
## MATLAB/Octave Implementation
Step 1. Load data
```matlab
data = readtable('global_mean_sea_level_1993-2024.csv');
valid_idx = ~isnan(data.GMSLWithGIA);
t = data.YearPlusFraction(valid_idx);
S = data.GMSLWithGIA(valid_idx);
t_centered = t - mean(t);
fprintf('Data points: %d\n', length(S));
```
Purpose: Load and clean GMSL data.
Step 2. Build Fourier design matrix
```matlab
n = length(t);
X_fourier = ones(n, 5);
X_fourier(:, 2) = cos(2*pi*1*t); % Annual cosine
X_fourier(:, 3) = sin(2*pi*1*t); % Annual sine
X_fourier(:, 4) = cos(2*pi*2*t); % Semi-annual cosine
X_fourier(:, 5) = sin(2*pi*2*t); % Semi-annual sine
```
Purpose: Create trigonometric basis functions.
Step 3. Solve for Fourier coefficients
```matlab
beta_fourier = (X_fourier' * X_fourier) \ (X_fourier' * S);
S_seasonal = X_fourier * beta_fourier;
fprintf('Fourier coefficients:\n');
disp(round(beta_fourier, 2));
```
Purpose: Extract seasonal component via [least squares](https://en.wikipedia.org/wiki/Least_squares).
Step 4. De-seasonalize and fit linear trend
```matlab
residuals = S - S_seasonal;
X_linear = [ones(size(t_centered)), t_centered];
beta_linear = (X_linear' * X_linear) \ (X_linear' * residuals);
trend_linear = X_linear * beta_linear;
fprintf('Sea level rise rate: %.3f mm/year\n', beta_linear(2));
```
Purpose: Estimate trend from de-seasonalized data.
Step 5. Fit quadratic trend
```matlab
X_quad = [ones(size(t_centered)), t_centered, t_centered.^2];
beta_quad = (X_quad' * X_quad) \ (X_quad' * residuals);
trend_quad = X_quad * beta_quad;
fprintf('Acceleration: %.4f mm/year²\n', 2*beta_quad(3));
```
Purpose: Test for acceleration.
Step 6. Calculate R² values
```matlab
SS_tot = sum((residuals - mean(residuals)).^2);
SS_res_linear = sum((residuals - trend_linear).^2);
SS_res_quad = sum((residuals - trend_quad).^2);
R2_linear = 1 - SS_res_linear/SS_tot;
R2_quad = 1 - SS_res_quad/SS_tot;
fprintf('R² linear: %.4f\n', R2_linear);
fprintf('R² quadratic: %.4f\n', R2_quad);
```
Purpose: Quantify model fit quality.
Step 7. Visualization
```matlab
figure('Position', [100, 100, 800, 600]);
% Original data
subplot(3,1,1);
plot(t, S, 'b-', 'LineWidth', 0.5);
ylabel('Sea Level (mm)'); title('Original GMSL Data');
grid on;
% Seasonal component
subplot(3,1,2);
plot(t, S_seasonal - mean(S_seasonal), 'g-', 'LineWidth', 1);
ylabel('Seasonal (mm)'); title('Seasonal Component');
grid on;
% De-seasonalized with trend
subplot(3,1,3);
plot(t, residuals, 'Color', [0.5 0.5 0.5], 'LineWidth', 0.5);
hold on;
plot(t, trend_quad, 'r-', 'LineWidth', 2);
xlabel('Year'); ylabel('Residual (mm)');
title('De-seasonalized Data with Quadratic Trend');
legend('Residuals', 'Trend'); grid on;
```
Purpose: Show complete decomposition.
Step 8. Final residual patterns
```matlab
final_residuals = residuals - trend_quad;
figure;
plot(t, final_residuals, 'k-', 'LineWidth', 0.5);
hold on; yline(0, 'r--');
xlabel('Year'); ylabel('Final Residual (mm)');
title('Residuals after Seasonal + Trend Removal');
grid on;
% Observe patterns corresponding to El Niño events
```
Purpose: Reveal climate oscillations in final residuals.
# Results Summary
## Fourier Coefficients (approximate)
- $a_0 \approx 14.5$ mm (mean level)
- $a_1 \approx 1.0$ mm (annual cosine)
- $b_1 \approx -4.8$ mm (annual sine)
- $a_2 \approx -1.6$ mm (semi-annual cosine)
- $b_2 \approx 0.4$ mm (semi-annual sine)
These coefficients represent the amplitude of [harmonic oscillations](https://en.wikipedia.org/wiki/Harmonic_oscillator) at specific frequencies.
## Trend Estimates
| Model | Approach | Rate (mm/yr) | Acceleration (mm/yr²) | R² |
|-------|----------|--------------|----------------------|-----|
| Linear | Direct | ~3.30 | — | ~0.964 |
| Linear | De-seasoned | ~3.30 | — | ~0.976 |
| Quadratic | Direct | ~3.30 | ~0.038 | ~0.973 |
| Quadratic | De-seasoned | ~3.30 | ~0.038 | ~0.985 |
**Key Finding**: Trend parameters are essentially identical whether or not you remove seasonality, but [R²](https://en.wikipedia.org/wiki/Coefficient_of_determination) improves significantly with [seasonal adjustment](https://en.wikipedia.org/wiki/Seasonal_adjustment).
# Physical Interpretation
## Seasonal Components
- **Annual cycle** (period = 1 year): Dominant seasonal pattern from [Earth's orbital mechanics](https://en.wikipedia.org/wiki/Earth%27s_orbit) and [solar irradiance](https://en.wikipedia.org/wiki/Solar_irradiance)
- **Semi-annual cycle** (period = 0.5 years): Secondary pattern from hemispheric asymmetry and [monsoon](https://en.wikipedia.org/wiki/Monsoon) cycles
## Trend Components
- **Linear rise**: ~3.3 mm/year average sea level rise (compare to [NOAA's estimate](https://www.climate.gov/news-features/understanding-climate/climate-change-global-sea-level))
- **Acceleration**: ~0.038 mm/year² suggests increasing rate
- **Physical drivers**: [Thermal expansion](https://en.wikipedia.org/wiki/Thermal_expansion) + [ice sheet melt](https://climate.nasa.gov/vital-signs/land-ice/)
## Residual Patterns
After removing seasonal and trend components, residuals show:
- **[El Niño](https://www.climate.gov/enso) events** (1997-98, 2015-16): Positive spikes ([ONI data](https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php))
- **[La Niña](https://en.wikipedia.org/wiki/La_Niña) events**: Negative deviations
- **Period**: Irregular 2-7 year cycles ([ENSO](https://en.wikipedia.org/wiki/El_Niño–Southern_Oscillation)) not captured by annual harmonics
# Reflection
## 1. Frequency Selection
**Question**: We used $k=1$ (annual) and $k=2$ (semi-annual). What would happen if you added $k=3$ (tri-annual)?
*Where to add this code: After Step 5 (solving for Fourier coefficients) in your implementation.*
### Python
```python
# Extended Fourier design matrix with k=3
X_fourier_extended = np.ones((len(t), 7))
X_fourier_extended[:, 1] = np.cos(2*np.pi*1*t) # k=1
X_fourier_extended[:, 2] = np.sin(2*np.pi*1*t)
X_fourier_extended[:, 3] = np.cos(2*np.pi*2*t) # k=2
X_fourier_extended[:, 4] = np.sin(2*np.pi*2*t)
X_fourier_extended[:, 5] = np.cos(2*np.pi*3*t) # k=3 (NEW)
X_fourier_extended[:, 6] = np.sin(2*np.pi*3*t)
# Solve and compare
beta_extended = np.linalg.solve(X_fourier_extended.T @ X_fourier_extended,
X_fourier_extended.T @ S)
print(f"k=3 coefficients: a3={beta_extended[5]:.3f}, b3={beta_extended[6]:.3f}")
# Expect very small coefficients - no 4-month cycle in sea level
```
### R
```r
# Extended design matrix with k=3
X_fourier_ext <- matrix(1, nrow=n, ncol=7)
X_fourier_ext[, 2] <- cos(2*pi*1*t) # k=1
X_fourier_ext[, 3] <- sin(2*pi*1*t)
X_fourier_ext[, 4] <- cos(2*pi*2*t) # k=2
X_fourier_ext[, 5] <- sin(2*pi*2*t)
X_fourier_ext[, 6] <- cos(2*pi*3*t) # k=3 (NEW)
X_fourier_ext[, 7] <- sin(2*pi*3*t)
# Solve and compare
beta_ext <- solve(t(X_fourier_ext) %*% X_fourier_ext) %*% (t(X_fourier_ext) %*% S)
cat("k=3 coefficients: a3=", round(beta_ext[6], 3),
" b3=", round(beta_ext[7], 3), "\n")
```
### MATLAB
```matlab
% Extended design matrix with k=3
X_fourier_ext = ones(n, 7);
X_fourier_ext(:, 2) = cos(2*pi*1*t); % k=1
X_fourier_ext(:, 3) = sin(2*pi*1*t);
X_fourier_ext(:, 4) = cos(2*pi*2*t); % k=2
X_fourier_ext(:, 5) = sin(2*pi*2*t);
X_fourier_ext(:, 6) = cos(2*pi*3*t); % k=3 (NEW)
X_fourier_ext(:, 7) = sin(2*pi*3*t);
% Solve and compare
beta_ext = (X_fourier_ext' * X_fourier_ext) \ (X_fourier_ext' * S);
fprintf('k=3 coefficients: a3=%.3f, b3=%.3f\n', beta_ext(6), beta_ext(7));
```
## 2. Orthogonality Check
**Question**: Verify that the Fourier basis functions are orthogonal over the data period.
*Where to add this code: After Step 4 (building the Fourier design matrix).*
### Python
```python
# Check orthogonality between different basis functions
ortho_12 = X_fourier[:, 1].T @ X_fourier[:, 2] # cos vs sin of same freq
ortho_13 = X_fourier[:, 1].T @ X_fourier[:, 3] # cos of different freqs
print(f"cos(2πt) · sin(2πt) = {ortho_12:.6f}")
print(f"cos(2πt) · cos(4πt) = {ortho_13:.6f}")
# Should be near zero (exact orthogonality only for integer periods)
```
### R
```r
# Check orthogonality
ortho_12 <- sum(X_fourier[, 2] * X_fourier[, 3]) # cos vs sin
ortho_13 <- sum(X_fourier[, 2] * X_fourier[, 4]) # different freqs
cat("cos(2πt) · sin(2πt) =", round(ortho_12, 6), "\n")
cat("cos(2πt) · cos(4πt) =", round(ortho_13, 6), "\n")
```
### MATLAB
```matlab
% Check orthogonality
ortho_12 = X_fourier(:, 2)' * X_fourier(:, 3); % cos vs sin
ortho_13 = X_fourier(:, 2)' * X_fourier(:, 4); % different freqs
fprintf('cos(2πt) · sin(2πt) = %.6f\n', ortho_12);
fprintf('cos(2πt) · cos(4πt) = %.6f\n', ortho_13);
```
## 3. Model Comparison
**Question**: The quadratic model has R² = 0.985 vs 0.976 for linear. Is the acceleration term statistically significant?
The [F-test](https://en.wikipedia.org/wiki/F-test) determines whether the quadratic model's better fit (R² =0.985 vs 0.976) is statistically significant or just due to having an extra parameter - it compares the improvement in fit against what we'd expect by random chance. A large [F-statistic](https://en.wikipedia.org/wiki/F-statistic) (>4) with small [p-value](https://en.wikipedia.org/wiki/P-value) (<0.05) means the acceleration term captures real structure in the data, confirming that sea level rise is genuinely accelerating rather than rising at a constant rate. The following code can be added to your working code to assess this.
*Where to add this code: After Step 10 (calculating R² values).*
### Python
```python
# F-test for nested models
from scipy import stats
n = len(residuals)
p_linear = 2 # parameters in linear model
p_quad = 3 # parameters in quadratic model
# Residual sum of squares
RSS_linear = np.sum((residuals - trend_linear)**2)
RSS_quad = np.sum((residuals - trend_quad)**2)
# F-statistic
F_stat = ((RSS_linear - RSS_quad)/(p_quad - p_linear)) / (RSS_quad/(n - p_quad))
# P-value
p_value = 1 - stats.f.cdf(F_stat, p_quad - p_linear, n - p_quad)
print(f"F-statistic: {F_stat:.2f}")
print(f"p-value: {p_value:.6f}")
if p_value < 0.05:
print("Acceleration is statistically significant at α=0.05")
else:
print("Acceleration is NOT significant at α=0.05")
```
### R
```r
# F-test for nested models
n <- length(residuals)
p_linear <- 2
p_quad <- 3
# F-statistic
F_stat <- ((SS_res_linear - SS_res_quad)/(p_quad - p_linear)) /
(SS_res_quad/(n - p_quad))
# P-value
p_value <- 1 - pf(F_stat, p_quad - p_linear, n - p_quad)
cat("F-statistic:", round(F_stat, 2), "\n")
cat("p-value:", format(p_value, scientific=FALSE, digits=6), "\n")
if (p_value < 0.05) {
cat("Acceleration is statistically significant at α=0.05\n")
} else {
cat("Acceleration is NOT significant at α=0.05\n")
}
```
### MATLAB
```matlab
% F-test for nested models
n = length(residuals);
p_linear = 2;
p_quad = 3;
% F-statistic
F_stat = ((SS_res_linear - SS_res_quad)/(p_quad - p_linear)) / ...
(SS_res_quad/(n - p_quad));
% P-value (using F-distribution CDF)
p_value = 1 - fcdf(F_stat, p_quad - p_linear, n - p_quad);
fprintf('F-statistic: %.2f\n', F_stat);
fprintf('p-value: %.6f\n', p_value);
if p_value < 0.05
fprintf('Acceleration is statistically significant at α=0.05\n');
else
fprintf('Acceleration is NOT significant at α=0.05\n');
end
```
## 4. Prediction
**Question**: Using the quadratic model, what sea level would you predict for 2050?
*Where to add this code: After Step 9 (fitting quadratic trend).* Note: See [IPCC projections](https://www.ipcc.ch/srocc/chapter/chapter-4-sea-level-rise-and-implications-for-low-lying-islands-coasts-and-communities/) for comparison.
### Python
```python
# Predict for 2050
t_2050 = 2050.5 # Mid-year 2050
t_2050_centered = t_2050 - t.mean()
prediction_2050 = (beta_quad[0] +
beta_quad[1] * t_2050_centered +
beta_quad[2] * t_2050_centered**2)
# Add back seasonal mean
prediction_2050 += beta_fourier[0]
print(f"Predicted sea level for 2050: {prediction_2050:.1f} mm")
print(f"Current (2024): ~{S[-1]:.1f} mm")
print(f"Projected rise: {prediction_2050 - S[-1]:.1f} mm")
# WARNING: Extrapolation assumes trends continue unchanged!
```
### R
```r
# Predict for 2050
t_2050 <- 2050.5
t_2050_centered <- t_2050 - mean(t)
prediction_2050 <- beta_quad[1] +
beta_quad[2] * t_2050_centered +
beta_quad[3] * t_2050_centered^2
# Add back seasonal mean
prediction_2050 <- prediction_2050 + beta_fourier[1]
cat("Predicted sea level for 2050:", round(prediction_2050, 1), "mm\n")
cat("Current (2024): ~", round(tail(S, 1), 1), "mm\n")
cat("Projected rise:", round(prediction_2050 - tail(S, 1), 1), "mm\n")
```
### MATLAB
```matlab
% Predict for 2050
t_2050 = 2050.5;
t_2050_centered = t_2050 - mean(t);
prediction_2050 = beta_quad(1) + ...
beta_quad(2) * t_2050_centered + ...
beta_quad(3) * t_2050_centered^2;
% Add back seasonal mean
prediction_2050 = prediction_2050 + beta_fourier(1);
fprintf('Predicted sea level for 2050: %.1f mm\n', prediction_2050);
fprintf('Current (2024): ~%.1f mm\n', S(end));
fprintf('Projected rise: %.1f mm\n', prediction_2050 - S(end));
```
## 5. Alternative Decomposition
**Question**: Instead of pre-specifying frequencies, how could you discover dominant frequencies?
*This would replace Steps 4-5 with frequency discovery first using [spectral analysis](https://en.wikipedia.org/wiki/Spectral_density).*
### Python
```python
# Use FFT to find dominant frequencies
from scipy.fft import fft, fftfreq
# Compute FFT of centered data
S_centered = S - S.mean()
yf = fft(S_centered)
xf = fftfreq(len(S), 1/12)[:len(S)//2] # 12 samples/year
# Find peaks in power spectrum
power = 2.0/len(S) * np.abs(yf[:len(S)//2])
peaks = np.where(power > np.percentile(power, 95))[0]
print("Dominant frequencies (cycles/year):", xf[peaks])
```
### R
```r
# FFT to find frequencies
S_centered <- S - mean(S)
ft <- fft(S_centered)
freq <- (0:(length(S)-1)) * 12 / length(S) # 12 samples/year
power <- Mod(ft)^2 / length(S)
# Find dominant frequencies (first half of spectrum)
half_len <- floor(length(S)/2)
dominant <- which(power[1:half_len] > quantile(power[1:half_len], 0.95))
cat("Dominant frequencies (cycles/year):", freq[dominant], "\n")
```
### MATLAB
```matlab
% FFT to find frequencies
S_centered = S - mean(S);
Y = fft(S_centered);
f = (0:length(S)-1) * 12 / length(S); % 12 samples/year
power = abs(Y).^2 / length(S);
% Find peaks
half_len = floor(length(S)/2);
threshold = prctile(power(1:half_len), 95);
dominant = find(power(1:half_len) > threshold);
fprintf('Dominant frequencies (cycles/year): ');
fprintf('%.2f ', f(dominant));
fprintf('\n');
```
# Key Takeaways
1. **Fourier regression = OLS**: Once frequencies are specified, it's just [linear regression](https://en.wikipedia.org/wiki/Linear_regression) with trigonometric features
2. **Design matrix structure**: Columns are cos/sin at different frequencies
3. **Seasonal removal improves R²**: But trend estimates remain stable
4. **Physical interpretation matters**: Frequencies should match known phenomena
5. **Residual analysis reveals hidden patterns**: Climate oscillations appear after removing known components
6. **Connection to course themes**: Combines [linear algebra](https://en.wikipedia.org/wiki/Linear_algebra) (normal equations), [eigenanalysis](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix) (from PCA), and [Fourier series](https://en.wikipedia.org/wiki/Fourier_series)
# Extensions for Exploration
1. **Add more frequencies**: What happens with $k = 3, 4, 5$?
2. **Variable amplitude**: Allow seasonal amplitude to change over time
3. **Spectral analysis**: Use [FFT](https://en.wikipedia.org/wiki/Fast_Fourier_transform) to find optimal frequencies
4. **Cross-validation**: Split data to test predictive power ([cross-validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics)))
5. **Confidence intervals**: [Bootstrap](https://en.wikipedia.org/wiki/Bootstrapping_(statistics)) to get uncertainty estimates
# Mathematical Connections
- **From [[MATH310S26-Day3-WorkdayMaterials (Two-by-two matrices, normal equations, and ordinary least squares)|Day 3]]**: [Normal equations](https://en.wikipedia.org/wiki/Normal_equation) $X^TX\beta = X^Ty$
- **From [[MATH310S26-Day6-Work|Day 6]]**: Variance decomposition (like [PCA](https://en.wikipedia.org/wiki/Principal_component_analysis) but in time)
- **From [[MATH310S26-Day7-Notes|Day 7]]**: Fourier basis functions
- **To future**: [FFT](https://en.wikipedia.org/wiki/Fast_Fourier_transform) for frequency discovery, [wavelet analysis](https://en.wikipedia.org/wiki/Wavelet_transform) for time-varying frequencies
# Code and Data Files
- **Data**: `global_mean_sea_level_1993-2024.csv` (1169 rows × multiple columns)
- **Key columns**: `YearPlusFraction` (time), `GMSLWithGIA` (sea level with [glacial isostatic adjustment](https://en.wikipedia.org/wiki/Post-glacial_rebound))
- **All implementations produce identical numerical results** (verified to 3 decimal places)
- **Original data source**: [AVISO+](https://www.aviso.altimetry.fr/en/data/products/ocean-indicators-products/mean-sea-level.html) satellite altimetry
---
*Remember: In real applications, always validate assumptions, check residuals for patterns, and consider physical interpretability of your models!*