From e7bdccfa7fe6b0ba85fc17590128066485ab6bc5 Mon Sep 17 00:00:00 2001
From: "Asthana, Shivanshi" <shivanshi.asthana@studium.uni-hamburg.de>
Date: Wed, 22 Nov 2023 21:17:09 +0000
Subject: [PATCH] READMEForCode.md
---
READMEForCode.md | 112 +++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 112 insertions(+)
create mode 100644 READMEForCode.md
diff --git a/READMEForCode.md b/READMEForCode.md
new file mode 100644
index 0000000..8f683b9
--- /dev/null
+++ b/READMEForCode.md
@@ -0,0 +1,112 @@
+import xarray as xr
+import numpy as np
+import matplotlib.pyplot as plt
+from statsmodels.tsa.arima.model import ARIMA
+from sklearn.metrics import mean_squared_error
+from scipy.stats import pearsonr
+
+def plot_arima_forecast(order, years, time_series, future_years, label, color, original_legend_added=False):
+ model = ARIMA(time_series, order=order)
+ model_fit = model.fit()
+ forecast = model_fit.forecast(steps=len(future_years))
+
+ # Add legend for original data only once
+ if not original_legend_added:
+ plt.plot(years, time_series, label='Original Data', marker='o', color='green')
+ original_legend_added = True
+
+ plt.axvline(x=2022, color='gray', linestyle='--') # Vertical line at the year 2022
+ plt.plot(future_years, forecast, label=label, linestyle='dashed', marker='o', color=color)
+
+ # Hindcasts for lead times 1, 2, 3, 4, and 5 years
+ for lead_time in [1, 2, 3, 4, 5]:
+ hindcast_years = range(1981 + lead_time, 2023)
+ hindcast = model_fit.predict(start=lead_time, end=len(years) - 1)
+ plt.plot(hindcast_years, hindcast, linestyle='dashed', marker='o', color=color)
+
+ # Calculate RMSE for each hindcast
+ rmse_values = []
+ for lead_time in [1, 2, 3, 4, 5]:
+ hindcast = model_fit.predict(start=lead_time, end=len(years) - 1)
+ rmse = np.sqrt(mean_squared_error(time_series[lead_time:], hindcast))
+ rmse_values.append(rmse)
+
+ # Calculate Pearson correlation coefficients
+ pearson_coefficients = []
+ for lead_time in [1, 2, 3, 4, 5]:
+ hindcast = model_fit.predict(start=lead_time, end=len(years) - 1)
+ pearson_coefficient, _ = pearsonr(time_series[lead_time:], hindcast)
+ pearson_coefficients.append(pearson_coefficient)
+
+ return rmse_values, pearson_coefficients, original_legend_added
+
+# Specify the range of years
+years = range(1981, 2023)
+
+# Initialize a list to store spatial averages
+spatial_averages = []
+
+# Loop through each year
+for year in years:
+ # Load the count file for the current year
+ count_data = xr.open_dataset(f'/home/u/u301871/counts/count_{year}.nc')
+ cluster_count = count_data['cluster_count'].values # Adjust variable name
+
+ # Calculate the spatial average
+ spatial_average = np.mean(cluster_count)
+
+ # Append the result to the list
+ spatial_averages.append(spatial_average)
+
+ # Close the count dataset to free up resources
+ count_data.close()
+
+# Create a time series
+time_series = np.floor(spatial_averages)
+
+# Plot ARIMA forecasts with different parameters
+future_years = range(2023, 2034)
+original_legend_added = False
+rmse_1, pearson_1, original_legend_added = plot_arima_forecast(order=(20, 1, 3), years=years, time_series=time_series, future_years=future_years, label='ARIMA Forecast (20, 1, 3)', color='blue', original_legend_added=original_legend_added)
+rmse_2, pearson_2, original_legend_added = plot_arima_forecast(order=(20, 1, 1), years=years, time_series=time_series, future_years=future_years, label='ARIMA Forecast (20, 1, 1)', color='black', original_legend_added=original_legend_added)
+rmse_3, pearson_3, original_legend_added = plot_arima_forecast(order=(20, 1, 5), years=years, time_series=time_series, future_years=future_years, label='ARIMA Forecast (20, 1, 5)', color='red', original_legend_added=original_legend_added)
+rmse_4, pearson_4, original_legend_added = plot_arima_forecast(order=(20, 1, 2), years=years, time_series=time_series, future_years=future_years, label='ARIMA Forecast (20, 1, 2)', color='purple', original_legend_added=original_legend_added)
+rmse_5, pearson_5, original_legend_added = plot_arima_forecast(order=(20, 1, 4), years=years, time_series=time_series, future_years=future_years, label='ARIMA Forecast (20, 1, 4)', color='yellow', original_legend_added=original_legend_added)
+
+# Add legend to the first plot
+plt.legend()
+
+# Save the first plot as timeseries.png
+plt.savefig('timeseries.png')
+
+# Plot RMSE values
+lead_times = [1, 2, 3, 4, 5]
+plt.figure()
+for i, rmse_values in enumerate([rmse_1, rmse_2, rmse_3, rmse_4, rmse_5]):
+ plt.plot(lead_times, rmse_values, label=f'ARIMA Forecast (20, 1, {i + 1})', marker='o')
+
+plt.xlabel('Lead Time (Years)')
+plt.ylabel('RMSE')
+plt.title('RMSE of ARIMA Forecasts for Different Lead Times')
+plt.legend()
+plt.grid(True)
+
+# Save the second plot as rmse.png
+plt.savefig('rmse.png')
+
+# Plot Pearson correlation coefficients
+plt.figure()
+for i, pearson_values in enumerate([pearson_1, pearson_2, pearson_3, pearson_4, pearson_5]):
+ plt.plot(lead_times, pearson_values, label=f'ARIMA Forecast (20, 1, {i + 1})', marker='o')
+
+plt.xlabel('Lead Time (Years)')
+plt.ylabel('Pearson Correlation Coefficient')
+plt.title('Pearson Correlation Coefficient of ARIMA Forecasts for Different Lead Times')
+plt.legend()
+plt.grid(True)
+
+# Save the third plot as pearson.png
+plt.savefig('pearson.png')
+
+# Show the plots
+plt.show()
--
GitLab