## Introduction Climate scientists track global sea levels to understand long-term trends. Financial analysts monitor stock prices to predict market movements. Medical researchers study patient data to identify treatment patterns. All of these professionals face the same fundamental challenge: how do we separate meaningful signals from random noise in time-varying data? Fourier regression provides a powerful mathematical lens for this problem. By decomposing complex signals into simple sine and cosine waves—called **Fourier modes**—we can identify dominant frequencies, remove unwanted noise, and build predictive models. But here's the key insight: Fourier analysis isn't just about frequencies—it's fundamentally a least squares problem in disguise. ## Project Description Starting with global sea level data spanning three decades, you'll discover how linear algebra transforms our understanding of temporal patterns. This project bridges time series analysis with the regression techniques you've already mastered. **1 Point - Foundation**: Implement Fourier regression using ordinary least squares to model sea level trends. Construct design matrices with trigonometric functions, solve normal equations, and interpret the physical meaning of discovered Fourier modes through systematic analysis of your coefficient estimates. **2 Points - Exploration**: Continue the analysis by examining the residuals to determine the annual sea-level increase and investigate whether the data suggests this increase is accelerating. Use the provided code implementations to systematically analyze trend components and report findings through embedded questions, quantitative results, and interpretive graphics. **Key Deliverable**: An analysis demonstrating how least squares principles underlie frequency domain analysis, with specific focus on extracting and interpreting trend components from sea level data to assess climate change indicators. ## Mathematical Background: Fourier Regression as Linear Algebra The key insight is that Fourier analysis can be formulated as a standard least squares problem. Instead of fitting $y = \beta_0 + \beta_1 x$, we fit a model of trigonometric functions—each representing a distinct **Fourier mode**. For more about least-squares regression see: [[MATH307Su25 - Least Squares Approximation for $(x,y)$ Scatter Data and Applications | Least Squares Approximation for $(x,y)$ Scatter Data and Applications]] ### The Fourier Regression Model Given time series data $y(t_1), y(t_2), \ldots, y(t_n)$, we seek to approximate the signal as: $y(t) = a_0 + \sum_{k=1}^{m} a_k \cos(2\pi f_k t) + b_k \sin(2\pi f_k t)$ where $f_k$ represents the $k$-th frequency component (or **mode**). This is a linear combination of known functions, making it a perfect candidate for least squares regression. ### Design Matrix Construction The design matrix for Fourier regression becomes: $X = \begin{bmatrix} 1 & \cos(2\pi f_1 t_1) & \sin(2\pi f_1 t_1) & \cos(2\pi f_2 t_1) & \sin(2\pi f_2 t_1) & \cdots \\ 1 & \cos(2\pi f_1 t_2) & \sin(2\pi f_1 t_2) & \cos(2\pi f_2 t_2) & \sin(2\pi f_2 t_2) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ 1 & \cos(2\pi f_1 t_n) & \sin(2\pi f_1 t_n) & \cos(2\pi f_2 t_n) & \sin(2\pi f_2 t_n) & \cdots \end{bmatrix}$ The parameter vector contains our Fourier coefficients: $\boldsymbol{\beta} = \begin{bmatrix} a_0 \\ a_1 \\ b_1 \\ a_2 \\ b_2 \\ \vdots \end{bmatrix}$ ### Worked Example: Sea Level Analysis Using the global mean sea level data from 1993-2024, we need to make intelligent choices about which frequencies to include in our model. #### Frequency Selection Strategy For climate data, we expect several natural periodicities: - **Annual cycle** ($f_1 = 1$ cycle per year): Seasonal warming/cooling affects thermal expansion - **Semi-annual cycle** ($f_2 = 2$ cycles per year): Solar heating patterns and atmospheric pressure variations - **Long-term trend**: Captured through polynomial terms or lower-frequency modes **Check Your Understanding**: Why do we choose $f_1 = 1$ cycle per year as our fundamental frequency? What physical processes in the ocean-atmosphere system would produce an annual signal in sea level measurements? #### Setting Up the Normal Equations For our two-mode model with annual and semi-annual components: $X = \begin{bmatrix} 1 & \cos(2\pi \cdot 1 \cdot t_1) & \sin(2\pi \cdot 1 \cdot t_1) & \cos(2\pi \cdot 2 \cdot t_1) & \sin(2\pi \cdot 2 \cdot t_1) \\ 1 & \cos(2\pi \cdot 1 \cdot t_2) & \sin(2\pi \cdot 1 \cdot t_2) & \cos(2\pi \cdot 2 \cdot t_2) & \sin(2\pi \cdot 2 \cdot t_2) \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & \cos(2\pi \cdot 1 \cdot t_n) & \sin(2\pi \cdot 1 \cdot t_n) & \cos(2\pi \cdot 2 \cdot t_n) & \sin(2\pi \cdot 2 \cdot t_n) \end{bmatrix}$ **Check Your Understanding**: In the sea level dataset, time is measured as "YearPlusFraction" (e.g., 1993.5 represents mid-1993). Why is this time format convenient for our frequency analysis? What would happen if we used time in different units? #### Computing the Solution From the normal equations $X^T X \hat{\boldsymbol{\beta}} = X^T \mathbf{y}$, we obtain: $X^T X = \begin{bmatrix} n & \sum\cos(2\pi t_i) & \sum\sin(2\pi t_i) & \sum\cos(4\pi t_i) & \sum\sin(4\pi t_i) \\ \sum\cos(2\pi t_i) & \sum\cos^2(2\pi t_i) & \sum\cos(2\pi t_i)\sin(2\pi t_i) & \cdots & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots \end{bmatrix}$ **Check Your Understanding**: When we have long time series with regularly spaced data, many of the off-diagonal terms in $X^T X$ become small. Why? What mathematical property of sine and cosine functions causes this? #### Interpretation of Results From our analysis of the sea level data, we find: - $a_0 = 14.55$ mm (mean sea level offset from reference) - $a_1 = 1.03$, $b_1 = -4.84$ (annual cycle coefficients) - $a_2 = -1.59$, $b_2 = 0.39$ (semi-annual cycle coefficients) **Mode Amplitude Calculation**: The amplitude of the annual mode is: $A_1 = \sqrt{a_1^2 + b_1^2} = \sqrt{1.03^2 + (-4.84)^2} \approx 4.95 \text{ mm}$ This tells us the sea level oscillates with about ±5mm annual variation, which aligns with known thermal expansion effects. **Check Your Understanding**: The semi-annual amplitude is $A_2 = \sqrt{(-1.59)^2 + (0.39)^2} \approx 1.64$ mm. Why is this smaller than the annual amplitude? What does this suggest about the relative importance of different climate forcing mechanisms? ### Residual Analysis for Trend Detection After removing the periodic components, the residuals $\mathbf{r} = \mathbf{y} - X\hat{\boldsymbol{\beta}}$ reveal the underlying trend. These residuals can then be analyzed using polynomial regression: $r(t) = c_0 + c_1 t + c_2 t^2 + \varepsilon(t)$ **Check Your Understanding**: Why do we analyze trends in the residuals rather than including polynomial terms directly in the original Fourier model? What does this two-step approach tell us about the separability of periodic and trend components? ## Implementation Guidelines ### Part 1: Foundation Implementation and Analysis #### Step 1: Data Preparation and Design Matrix Construction **MATLAB Implementation:** ```matlab % Load and prepare sea level data data = readtable('global_mean_sea_level_19932024.csv'); time = data.YearPlusFraction; % Time in years sea_level = data.SmoothedGMSLWithGIA; % Sea level in mm % Remove any NaN values and display basic statistics valid_idx = ~isnan(sea_level); time = time(valid_idx); sea_level = sea_level(valid_idx); fprintf('Dataset spans %.2f to %.2f (%.1f years)\n', ... min(time), max(time), max(time) - min(time)); fprintf('Sea level range: %.1f to %.1f mm\n', min(sea_level), max(sea_level)); % Construct design matrix for 2-mode model n = length(time); f1 = 1; f2 = 2; % Annual and semi-annual frequencies (cycles per year) X = [ones(n,1), ... % Constant term cos(2*pi*f1*time), sin(2*pi*f1*time), ... % Annual mode cos(2*pi*f2*time), sin(2*pi*f2*time)]; % Semi-annual mode fprintf('Design matrix dimensions: %d x %d\n', size(X)); ``` **Python Implementation:** ```python import numpy as np import pandas as pd import matplotlib.pyplot as plt # Load and prepare sea level data data = pd.read_csv('global_mean_sea_level_19932024.csv') time = data['YearPlusFraction'].values # Time in years sea_level = data['SmoothedGMSLWithGIA'].values # Sea level in mm # Remove any NaN values and display basic statistics valid_mask = ~np.isnan(sea_level) time = time[valid_mask] sea_level = sea_level[valid_mask] print(f'Dataset spans {time.min():.2f} to {time.max():.2f} ({time.max() - time.min():.1f} years)') print(f'Sea level range: {sea_level.min():.1f} to {sea_level.max():.1f} mm') # Construct design matrix for 2-mode model n = len(time) f1, f2 = 1, 2 # Annual and semi-annual frequencies (cycles per year) X = np.column_stack([ np.ones(n), # Constant term np.cos(2*np.pi*f1*time), np.sin(2*np.pi*f1*time), # Annual mode np.cos(2*np.pi*f2*time), np.sin(2*np.pi*f2*time) # Semi-annual mode ]) print(f'Design matrix dimensions: {X.shape}') ``` #### Step 2: Solve Normal Equations and Analyze Coefficients **MATLAB:** ```matlab % Solve normal equations using matrix methods XtX = X' * X; % Compute X^T X Xty = X' * sea_level; % Compute X^T y % Check condition number before solving cond_num = cond(XtX); fprintf('Condition number of X^T X: %.2e\n', cond_num); % Solve for coefficients beta_hat = XtX \ Xty; % Solve normal equations y_fitted = X * beta_hat; % Compute fitted values % Extract and interpret coefficients a0 = beta_hat(1); % Mean offset a1 = beta_hat(2); b1 = beta_hat(3); % Annual mode coefficients a2 = beta_hat(4); b2 = beta_hat(5); % Semi-annual mode coefficients fprintf('\nFourier Regression Results:\n'); fprintf('Mean sea level offset: %.2f mm\n', a0); fprintf('Annual mode: a1 = %.3f, b1 = %.3f\n', a1, b1); fprintf('Semi-annual mode: a2 = %.3f, b2 = %.3f\n', a2, b2); % Calculate mode amplitudes and phases A1 = sqrt(a1^2 + b1^2); % Annual amplitude A2 = sqrt(a2^2 + b2^2); % Semi-annual amplitude phi1 = atan2(b1, a1); % Annual phase phi2 = atan2(b2, a2); % Semi-annual phase fprintf('\nMode Analysis:\n'); fprintf('Annual amplitude: %.2f mm, phase: %.2f radians\n', A1, phi1); fprintf('Semi-annual amplitude: %.2f mm, phase: %.2f radians\n', A2, phi2); ``` **Python:** ```python # Solve normal equations using matrix methods XtX = X.T @ X # Compute X^T X Xty = X.T @ sea_level # Compute X^T y # Check condition number before solving cond_num = np.linalg.cond(XtX) print(f'Condition number of X^T X: {cond_num:.2e}') # Solve for coefficients beta_hat = np.linalg.solve(XtX, Xty) # Solve normal equations y_fitted = X @ beta_hat # Compute fitted values # Extract and interpret coefficients a0 = beta_hat[0] # Mean offset a1, b1 = beta_hat[1], beta_hat[2] # Annual mode coefficients a2, b2 = beta_hat[3], beta_hat[4] # Semi-annual mode coefficients print('\nFourier Regression Results:') print(f'Mean sea level offset: {a0:.2f} mm') print(f'Annual mode: a1 = {a1:.3f}, b1 = {b1:.3f}') print(f'Semi-annual mode: a2 = {a2:.3f}, b2 = {b2:.3f}') # Calculate mode amplitudes and phases A1 = np.sqrt(a1**2 + b1**2) # Annual amplitude A2 = np.sqrt(a2**2 + b2**2) # Semi-annual amplitude phi1 = np.arctan2(b1, a1) # Annual phase phi2 = np.arctan2(b2, a2) # Semi-annual phase print('\nMode Analysis:') print(f'Annual amplitude: {A1:.2f} mm, phase: {phi1:.2f} radians') print(f'Semi-annual amplitude: {A2:.2f} mm, phase: {phi2:.2f} radians') ``` #### Step 3: Model Assessment and Visualization **MATLAB:** ```matlab % Calculate model fit statistics residuals = sea_level - y_fitted; RSS = sum(residuals.^2); TSS = sum((sea_level - mean(sea_level)).^2); R_squared = 1 - RSS/TSS; fprintf('\nModel Fit Assessment:\n'); fprintf('R-squared: %.4f\n', R_squared); fprintf('Root mean square error: %.2f mm\n', sqrt(RSS/n)); % Create comprehensive visualization figure('Position', [100, 100, 1400, 800]); % Time series plot subplot(2,2,1); plot(time, sea_level, 'b-', 'LineWidth', 1); hold on; plot(time, y_fitted, 'r-', 'LineWidth', 2); xlabel('Year'); ylabel('Sea Level (mm)'); title(sprintf('Sea Level Data and Fourier Fit (R² = %.3f)', R_squared)); legend('Observed', 'Fourier Model', 'Location', 'northwest'); % Residuals analysis subplot(2,2,2); plot(time, residuals, 'k-', 'LineWidth', 1); xlabel('Year'); ylabel('Residuals (mm)'); title('Residuals: Trend Component'); grid on; % Mode decomposition subplot(2,2,3); annual_component = a1*cos(2*pi*f1*time) + b1*sin(2*pi*f1*time); semiannual_component = a2*cos(2*pi*f2*time) + b2*sin(2*pi*f2*time); plot(time, annual_component, 'g-', 'LineWidth', 2); hold on; plot(time, semiannual_component, 'm-', 'LineWidth', 2); xlabel('Year'); ylabel('Amplitude (mm)'); title('Individual Fourier Modes'); legend('Annual', 'Semi-annual', 'Location', 'best'); % Frequency content visualization subplot(2,2,4); frequencies = [f1, f2]; amplitudes = [A1, A2]; bar(frequencies, amplitudes); xlabel('Frequency (cycles/year)'); ylabel('Amplitude (mm)'); title('Mode Amplitudes'); ``` **Python:** ```python # Calculate model fit statistics residuals = sea_level - y_fitted RSS = np.sum(residuals**2) TSS = np.sum((sea_level - np.mean(sea_level))**2) R_squared = 1 - RSS/TSS print('\nModel Fit Assessment:') print(f'R-squared: {R_squared:.4f}') print(f'Root mean square error: {np.sqrt(RSS/n):.2f} mm') # Create comprehensive visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Time series plot axes[0,0].plot(time, sea_level, 'b-', linewidth=1, label='Observed') axes[0,0].plot(time, y_fitted, 'r-', linewidth=2, label='Fourier Model') axes[0,0].set_xlabel('Year') axes[0,0].set_ylabel('Sea Level (mm)') axes[0,0].set_title(f'Sea Level Data and Fourier Fit (R² = {R_squared:.3f})') axes[0,0].legend(loc='upper left') # Residuals analysis axes[0,1].plot(time, residuals, 'k-', linewidth=1) axes[0,1].set_xlabel('Year') axes[0,1].set_ylabel('Residuals (mm)') axes[0,1].set_title('Residuals: Trend Component') axes[0,1].grid(True) # Mode decomposition annual_component = a1*np.cos(2*np.pi*f1*time) + b1*np.sin(2*np.pi*f1*time) semiannual_component = a2*np.cos(2*np.pi*f2*time) + b2*np.sin(2*np.pi*f2*time) axes[1,0].plot(time, annual_component, 'g-', linewidth=2, label='Annual') axes[1,0].plot(time, semiannual_component, 'm-', linewidth=2, label='Semi-annual') axes[1,0].set_xlabel('Year') axes[1,0].set_ylabel('Amplitude (mm)') axes[1,0].set_title('Individual Fourier Modes') axes[1,0].legend() # Frequency content visualization frequencies = [f1, f2] amplitudes = [A1, A2] axes[1,1].bar(frequencies, amplitudes) axes[1,1].set_xlabel('Frequency (cycles/year)') axes[1,1].set_ylabel('Amplitude (mm)') axes[1,1].set_title('Mode Amplitudes') plt.tight_layout() plt.show() ``` **Check Your Understanding**: Look at your residuals plot. Do you see a clear trend? What type of polynomial would best capture this pattern? How does this relate to discussions of sea level acceleration due to climate change? ### Part 2: Residual Analysis for Climate Trend Assessment **Primary Focus: Sea Level Rise Analysis** After removing periodic components, systematically analyze residuals to extract climate signals: **MATLAB:** ```matlab % Fit polynomial trends to residuals from Fourier model % Linear trend: measures constant rate of sea level rise X_linear = [ones(n,1), time]; beta_linear = (X_linear'*X_linear) \ (X_linear'*residuals); % Solves normal equations linear_fit = X_linear * beta_linear; RSS_linear = sum((residuals - linear_fit).^2); R2_linear = 1 - RSS_linear/sum((residuals - mean(residuals)).^2); % Quadratic trend: tests for acceleration in sea level rise X_quad = [ones(n,1), time, time.^2]; beta_quad = (X_quad'*X_quad) \ (X_quad'*residuals); % Solves normal equations quad_fit = X_quad * beta_quad; RSS_quad = sum((residuals - quad_fit).^2); R2_quad = 1 - RSS_quad/sum((residuals - mean(residuals)).^2); % Extract climate indicators annual_rise_rate = beta_linear(2); % mm per year acceleration = 2*beta_quad(3); % mm per year² fprintf('\n=== CLIMATE TREND ANALYSIS ===\n'); fprintf('Linear trend model:\n'); fprintf(' Annual sea level rise: %.2f ± %.2f mm/year\n', annual_rise_rate, std_error_linear); fprintf(' R²: %.4f\n', R2_linear); fprintf('\nQuadratic trend model:\n'); fprintf(' Sea level acceleration: %.3f mm/year²\n', acceleration); fprintf(' R²: %.4f\n', R2_quad); % Statistical significance test for acceleration if R2_quad > R2_linear + 0.01 % Threshold for meaningful improvement fprintf('\nAcceleration appears statistically significant\n'); fprintf('Quadratic model explains %.1f%% more variance\n', 100*(R2_quad-R2_linear)); else fprintf('\nNo clear evidence of acceleration in this dataset\n'); end ``` **Python:** ```python # Fit polynomial trends to residuals from Fourier model # Linear trend: measures constant rate of sea level rise X_linear = np.column_stack([np.ones(n), time]) beta_linear = np.linalg.solve(X_linear.T @ X_linear, X_linear.T @ residuals) linear_fit = X_linear @ beta_linear RSS_linear = np.sum((residuals - linear_fit)**2) R2_linear = 1 - RSS_linear/np.sum((residuals - np.mean(residuals))**2) # Quadratic trend: tests for acceleration in sea level rise X_quad = np.column_stack([np.ones(n), time, time**2]) beta_quad = np.linalg.solve(X_quad.T @ X_quad, X_quad.T @ residuals) quad_fit = X_quad @ beta_quad RSS_quad = np.sum((residuals - quad_fit)**2) R2_quad = 1 - RSS_quad/np.sum((residuals - np.mean(residuals))**2) # Extract climate indicators annual_rise_rate = beta_linear[1] # mm per year acceleration = 2*beta_quad[2] # mm per year² print('\n=== CLIMATE TREND ANALYSIS ===') print('Linear trend model:') print(f' Annual sea level rise: {annual_rise_rate:.2f} mm/year') print(f' R²: {R2_linear:.4f}') print('\nQuadratic trend model:') print(f' Sea level acceleration: {acceleration:.3f} mm/year²') print(f' R²: {R2_quad:.4f}') # Statistical significance test for acceleration if R2_quad > R2_linear + 0.01: # Threshold for meaningful improvement print('\nAcceleration appears statistically significant') print(f'Quadratic model explains {100*(R2_quad-R2_linear):.1f}% more variance') else: print('\nNo clear evidence of acceleration in this dataset') ``` **Check Your Understanding**: Based on your results, what is the annual rate of sea level rise? How does this compare to commonly cited values of ~3.3 mm/year? Is there evidence for acceleration in the rate of rise? **Visualization of Climate Trends** **MATLAB:** ```matlab % Create comprehensive trend analysis visualization figure('Position', [100, 100, 1200, 800]); % Subplot 1: Residuals with trend fits subplot(2,2,1); plot(time, residuals, 'k.', 'MarkerSize', 3); hold on; plot(time, linear_fit, 'b-', 'LineWidth', 2); plot(time, quad_fit, 'r-', 'LineWidth', 2); xlabel('Year'); ylabel('Residuals (mm)'); title('Sea Level Trends After Removing Seasonal Cycles'); legend('Residuals', 'Linear Trend', 'Quadratic Trend', 'Location', 'northwest'); grid on; % Subplot 2: Comparison of trend models subplot(2,2,2); plot(time, residuals, 'k-', 'LineWidth', 1, 'Color', [0.7 0.7 0.7]); hold on; plot(time, linear_fit, 'b-', 'LineWidth', 2); plot(time, quad_fit, 'r--', 'LineWidth', 2); xlabel('Year'); ylabel('Sea Level Trend (mm)'); title(sprintf('Linear vs Quadratic: %.2f vs %.2f mm/yr²', 0, acceleration)); legend('Detrended Data', 'Linear', 'Quadratic', 'Location', 'northwest'); % Subplot 3: Rate of change over time (for quadratic model) if abs(acceleration) > 0.001 % Only if acceleration is meaningful subplot(2,2,3); instantaneous_rate = beta_quad(2) + 2*beta_quad(3)*time; plot(time, instantaneous_rate, 'g-', 'LineWidth', 2); xlabel('Year'); ylabel('Instantaneous Rate (mm/year)'); title('Sea Level Rise Rate Over Time'); grid on; end % Subplot 4: Model comparison statistics subplot(2,2,4); models = {'Linear', 'Quadratic'}; r_squared_comparison = [R2_linear, R2_quad]; bar(r_squared_comparison); set(gca, 'XTickLabel', models); ylabel('R²'); title('Model Performance Comparison'); for i = 1:length(r_squared_comparison) text(i, r_squared_comparison(i)+0.005, sprintf('%.3f', r_squared_comparison(i)), ... 'HorizontalAlignment', 'center'); end ``` **Python:** ```python # Create comprehensive trend analysis visualization fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # Subplot 1: Residuals with trend fits axes[0,0].plot(time, residuals, 'k.', markersize=3, label='Residuals') axes[0,0].plot(time, linear_fit, 'b-', linewidth=2, label='Linear Trend') axes[0,0].plot(time, quad_fit, 'r-', linewidth=2, label='Quadratic Trend') axes[0,0].set_xlabel('Year') axes[0,0].set_ylabel('Residuals (mm)') axes[0,0].set_title('Sea Level Trends After Removing Seasonal Cycles') axes[0,0].legend(loc='upper left') axes[0,0].grid(True) # Subplot 2: Comparison of trend models axes[0,1].plot(time, residuals, 'k-', linewidth=1, color='lightgray', label='Detrended Data') axes[0,1].plot(time, linear_fit, 'b-', linewidth=2, label='Linear') axes[0,1].plot(time, quad_fit, 'r--', linewidth=2, label='Quadratic') axes[0,1].set_xlabel('Year') axes[0,1].set_ylabel('Sea Level Trend (mm)') axes[0,1].set_title(f'Linear vs Quadratic: 0 vs {acceleration:.2f} mm/yr²') axes[0,1].legend(loc='upper left') # Subplot 3: Rate of change over time (for quadratic model) if abs(acceleration) > 0.001: # Only if acceleration is meaningful instantaneous_rate = beta_quad[1] + 2*beta_quad[2]*time axes[1,0].plot(time, instantaneous_rate, 'g-', linewidth=2) axes[1,0].set_xlabel('Year') axes[1,0].set_ylabel('Instantaneous Rate (mm/year)') axes[1,0].set_title('Sea Level Rise Rate Over Time') axes[1,0].grid(True) # Subplot 4: Model comparison statistics models = ['Linear', 'Quadratic'] r_squared_comparison = [R2_linear, R2_quad] bars = axes[1,1].bar(models, r_squared_comparison) axes[1,1].set_ylabel('R²') axes[1,1].set_title('Model Performance Comparison') for i, v in enumerate(r_squared_comparison): axes[1,1].text(i, v+0.005, f'{v:.3f}', ha='center') plt.tight_layout() plt.show() ``` ## Check Your Understanding If you observe a positive acceleration coefficient, what does this suggest about future sea level projections? How might this information influence coastal planning and climate policy decisions? --- ```python # Calculate AIC and BIC values for Fourier models bic_values = np.zeros(len(max_modes)) aic_values = np.zeros(len(max_modes)) for i, num_modes in enumerate(max_modes): # Construct design matrix with specified number of modes X_modes = np.ones((n, 1)) # Start with constant term for k in range(1, num_modes + 1): freq = k # Use integer frequencies: 1, 2, 3, ... cycles/year cos_col = np.cos(2*np.pi*freq*time).reshape(-1, 1) sin_col = np.sin(2*np.pi*freq*time).reshape(-1, 1) X_modes = np.hstack([X_modes, cos_col, sin_col]) # Fit model using normal equations beta_modes = np.linalg.solve(X_modes.T @ X_modes, X_modes.T @ sea_level) fitted_modes = X_modes @ beta_modes # Calculate information criteria num_params = X_modes.shape[1] RSS = np.sum((sea_level - fitted_modes)**2) # AIC = n*log(RSS/n) + 2*p aic_values[i] = n*np.log(RSS/n) + 2*num_params # BIC = n*log(RSS/n) + log(n)*p bic_values[i] = n*np.log(RSS/n) + np.log(n)*num_params print(f'Modes: {num_modes}, Parameters: {num_params}, AIC: {aic_values[i]:.1f}, BIC: {bic_values[i]:.1f}') # Find optimal number of modes best_aic = np.argmin(aic_values) best_bic = np.argmin(bic_values) print('\nOptimal model selection:') print(f'Best AIC: {max_modes[best_aic]} modes') print(f'Best BIC: {max_modes[best_bic]} modes') ``` --- #### Student Tasks For the residual analysis: 1. **Implement** the trend analysis using matrix methods to extract linear and quadratic components 2. **Quantify** the annual sea level rise rate and acceleration coefficient from your regression results 3. **Assess** statistical significance by comparing $R^2$ values and model performance 4. **Interpret** the climate implications of your findings in the context of global sea level monitoring --- ## Analysis Framework ### Mathematical Insights - How does the **[condition number](Condition%20Number.md)** of the design matrix change as you add polynomial terms to the trend analysis? - What is the relationship between the Fourier model residuals and the underlying physics of sea level change? ### Climate Science Applications - How do your calculated rise rates compare to published IPCC estimates? - What does the presence or absence of acceleration tell us about climate change impacts? ### Statistical Considerations - How does the $R^2$ improvement from linear to quadratic models help assess the evidence for acceleration? - What does the residual structure tell us about unmodeled climate processes? **Check Your Understanding:** If your analysis shows significant acceleration, what physical processes could explain this trend? Consider thermal expansion, ice sheet melting, and gravitational effects. --- ## Real-World Context Understanding how to decompose signals into trend and periodic components has practical implications across many fields: ### Climate Science Applications - **[Sea Level Monitoring](https://climate.nasa.gov/evidence/):** Separating seasonal cycles from long-term trends to assess climate change impacts - **[Temperature Analysis](https://www.google.com/search?q=climate+temperature+trend+analysis+seasonal+adjustment):** Distinguishing natural variability from anthropogenic warming signals - **[Ice Sheet Dynamics](https://en.wikipedia.org/wiki/Ice_sheet_mass_balance):** Analyzing seasonal accumulation/ablation cycles versus long-term mass changes ### Oceanographic Applications - **[Tidal Analysis](https://en.wikipedia.org/wiki/Tidal_analysis):** Separating multiple harmonic components in sea level observations - **[Ocean Current Monitoring](https://www.google.com/search?q=oceanographic+time+series+analysis+current+measurements):** Identifying periodic flow patterns and long-term circulation changes - **[Marine Ecosystem Studies](https://en.wikipedia.org/wiki/Marine_ecosystem):** Understanding seasonal breeding cycles and population dynamics ### Policy and Planning Applications - **[Coastal Management](https://www.google.com/search?q=coastal+planning+sea+level+rise+projections):** Using trend analysis to inform infrastructure planning and zoning decisions - **[Climate Adaptation](https://en.wikipedia.org/wiki/Climate_change_adaptation):** Translating mathematical trends into actionable policy recommendations - **[Risk Assessment](https://www.google.com/search?q=coastal+flood+risk+sea+level+projections):** Quantifying future flood risks based on acceleration patterns ### Technical Challenges The intersection of time series analysis and climate science presents several considerations: - **[Measurement Uncertainty](https://en.wikipedia.org/wiki/Satellite_altimetry):** Understanding how instrument precision affects trend detection - **[Regional Variation](https://www.google.com/search?q=regional+sea+level+variation+global+mean):** Recognizing that global mean trends may not reflect local changes - **[Natural Variability](https://en.wikipedia.org/wiki/El_Ni%C3%B1o%E2%80%93Southern_Oscillation):** Distinguishing long-term climate signals from multi-year natural cycles ### Research Impact The mathematical techniques you're applying enable critical assessments of climate change. The ability to detect acceleration in sea level rise provides evidence for policy discussions about emission reductions and coastal adaptation. These methods help scientists distinguish between natural climate variability and anthropogenic changes, forming the foundation for climate projections used in international assessments. --- ## Deliverable Checklist - [ ] **Code Implementation** - [ ] Fourier regression using normal equations with proper trigonometric design matrix construction - [ ] Residual analysis using polynomial trend fitting with linear and quadratic models - [ ] Statistical assessment of model performance and trend significance - [ ] Clear documentation of coefficient interpretation and climate indicators - [ ] **Analysis Report** - [ ] Verification that Fourier model successfully removes seasonal cycles from sea level data - [ ] Quantification of annual sea level rise rate with uncertainty assessment - [ ] Investigation of acceleration evidence using quadratic trend analysis - [ ] Discussion of condition number effects and numerical stability in trend fitting - [ ] Physical interpretation of trend findings in context of climate change science - [ ] **Communication** - [ ] One-page synthesis connecting mathematical trend analysis to climate policy implications - [ ] Five to seven minute video walkthrough demonstrating progression from raw data to climate indicators - [ ] Clear explanation of how matrix-based frequency analysis enables separation of natural cycles from anthropogenic trends in climate data ---