import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import spearmanr
from numpy import random
# ------------------------------
# 1. Generate realistic rainfall data with more dry days
# ------------------------------
np.random.seed(42)
# Generate dry and wet spell lengths using a geometric distribution
dry_spell_lengths = np.random.geometric(p=0.9, size=50) # More frequent dry spells
wet_spell_lengths = np.random.geometric(p=0.3, size=10) # Less frequent wet spells
# Construct rainfall time series
rainfall_data = []
for dry, wet in zip(dry_spell_lengths, wet_spell_lengths):
rainfall_data.extend([0] * dry) # Add dry days
rainfall_data.extend(np.random.randint(0, 100, size=wet)) # Add wet days
(moderate rainfall)
rainfall_data = np.array(rainfall_data)
print(rainfall_data, len(rainfall_data))
# ------------------------------
# 2. Compute transition probabilities
# ------------------------------
n = len(rainfall_data) - 1
dry_prev = rainfall_data[:-1] == 0
wet_prev = rainfall_data[:-1] > 0
dry_curr = rainfall_data[1:] == 0
wet_curr = rainfall_data[1:] > 0
P_00 = np.sum(dry_prev & dry_curr) / n
P_01 = np.sum(dry_prev & wet_curr) / n
P_10 = np.sum(wet_prev & dry_curr) / n
P_11 = np.sum(wet_prev & wet_curr) / n
print(f"P_00 (Dry-Dry): {P_00:.3f}")
print(f"P_01 (Dry-Wet): {P_01:.3f}")
print(f"P_10 (Wet-Dry): {P_10:.3f}")
print(f"P_11 (Wet-Wet): {P_11:.3f}")
# ------------------------------
# 3. Compute empirical CDFs
# ------------------------------
def empirical_cdf(data, x):
return np.count_nonzero(data <= x) / len(data)
H_x = lambda x: empirical_cdf(rainfall_data[:-1], x)
H_y = lambda y: empirical_cdf(rainfall_data[1:], y)
# Conditional CDFs for wet days
wet_indices = (rainfall_data[:-1] > 0) & (rainfall_data[1:] > 0)
wet_x_data = rainfall_data[:-1][wet_indices]
wet_y_data = rainfall_data[1:][wet_indices]
# Estimate copula parameter using Spearman’s correlation
if len(wet_x_data) > 1:
rho, _ = spearmanr(wet_x_data, wet_y_data)
theta = 1 / (1 - rho) # Gumbel copula parameter
else:
theta = 1 # Default to independence if insufficient data
# ------------------------------
# 4. Gumbel Copula Functions
# ------------------------------
# Gumbel copula CDF
def gumbel_copula_cdf(u, v, theta):
if theta == 1:
return u * v # Independent case
inner_term = ((-np.log(u))**theta + (-np.log(v))**theta)**(1/theta)
return np.exp(-inner_term)
# Partial derivative of Gumbel copula w.r.t. u
def gumbel_copula_partial_u(u, v, theta):
if theta == 1:
return v # Independent case
inner_term = (-np.log(u))**theta + (-np.log(v))**theta
C_uv = gumbel_copula_cdf(u, v, theta)
return C_uv * ((-np.log(u))**(theta-1)) / (u * inner_term**(1 - 1/theta))
# Conditional CDF using Gumbel copula
def O_Y_given_X_x_and_X_positive(y, x):
u, v = H_x(x), H_y(y)
C_uv_partial = gumbel_copula_partial_u(u, v, theta)
h_x = np.count_nonzero(rainfall_data[:-1] == x) / len(rainfall_data[:-1]) if
len(rainfall_data[:-1]) > 0 else 0
f_y = np.count_nonzero(rainfall_data[1:] == y) / len(rainfall_data[1:]) if
len(rainfall_data[1:]) > 0 else 0
denom = P_10 * h_x + P_11 * f_y
return (P_10 * h_x + P_11 * f_y * C_uv_partial) / denom if denom > 0 else
H_y(y)
# Compute Conditional CDF for X = 0
def O_Y_given_X_0(y):
denom = P_00 + P_01
return (P_00 + P_01 * H_y(y)) / denom if denom > 0 else H_y(y)
# ------------------------------
# 5. Generate Plots for Conditional CDF
# ------------------------------
y_values = np.linspace(0, max(rainfall_data), 100)
cdf_0 = np.array([O_Y_given_X_0(y) for y in y_values])
cdf_x_positive = np.array([O_Y_given_X_x_and_X_positive(y, 3) for y in y_values])
plt.figure(figsize=(10, 6))
plt.plot(y_values, cdf_0, label="O_Y_given_X=0", linestyle="--", linewidth=2)
plt.plot(y_values, cdf_x_positive, label="O_Y_given_X=x and X>0", linewidth=2)
plt.scatter(rainfall_data[1:], [O_Y_given_X_x_and_X_positive(y, 3) for y in
rainfall_data[1:]],
color='red', alpha=0.5, label="Observed Rainfall")
plt.xlabel('Rainfall (Y)')
plt.ylabel('Conditional CDF')
plt.legend()
plt.title('Conditional CDF of Y given X')
plt.grid()
plt.show()
# ------------------------------
# 6. Simulate rainfall using the Markov Chain with Copula
# ------------------------------
simulated_rainfall = np.zeros(n, dtype=float)
# For the simulation, use the observed first day as the starting point
simulated_rainfall[0] = rainfall_data[0]
# For wet days, sample intensities from the observed wet days
wet_intensities = rainfall_data[rainfall_data > 0]
if len(wet_intensities) == 0:
raise ValueError("No wet days found in observed data; cannot sample
intensities.")
# Simulate each day based on the previous day's state:
for i in range(1, n):
prev = simulated_rainfall[i-1]
r = np.random.rand() # uniform random number in [0,1)
if prev == 0:
# Previous day was dry: use transitions from a dry day.
if r < P_00:
simulated_rainfall[i] = 0
else:
simulated_rainfall[i] = np.random.choice(wet_intensities)
else:
# Previous day was wet: use transitions from a wet day.
if r < P_10:
simulated_rainfall[i] = 0
else:
simulated_rainfall[i] = np.random.choice(wet_intensities)
# Plot simulated vs observed rainfall
plt.figure(figsize=(10, 6))
plt.plot(rainfall_data, label="Observed Rainfall", marker="o", linestyle="-",
alpha=0.7)
plt.plot(simulated_rainfall, label="Simulated Rainfall", marker="x",
linestyle="--", alpha=0.7)
plt.xlabel("Time (days)")
plt.ylabel("Rainfall (mm)")
plt.title("Observed vs Simulated Rainfall")
plt.legend()
plt.grid(True)
plt.show()
