Generate multivariate distributions of lognormal and normal distribution in python

Question

I need to generate random numbers from 3 correlated distributions. First two of them are lognormal and the final one is normal, i.e. for X, Y and Z I need X and Y to be lognormal, and Z is normal. The correlation matrix is given:

1, -0.8, 0.5
-0.8, 1,   -0.5
0.5,-0.5,  1

For example, the std of X and Y and Z are 0.5, 0.3 and 20. The mean of X and Y can be 1, and 0 for Z. How can I do that?

I know I can use https://stackoverflow.com/questions/16024677/generate-correlated-data-in-python-3-3 to generate 3 normally distributed. If all 3 are lognormal distributions, just simply use np.exp. But how can I apply it to my case?

Now my solution is: we can generate multivariate normal distribution of log(X), log(Y) and Z since these 3 variables are normal. However, we need to carefully calculate the elements of the covariance matrix and the vector of mean to ensure the requirements of the mean and std of X and Y and Z.

For mean vector, it is simply [-np.log(1 + std_X ** 2)/2, -np.log(1 + std_Y ** 2)/2, 0] since EX=EY=1. For covariance matrix, based on my calculation, the covariance between \ln X and \ln Y comes from $cov(\ln X, \ln Y)=\ln (1 + \frac{cov(X,Y)}{EXEY}) $. The covariance between \ln X and Z comes from https://stats.stackexchange.com/a/470146/303835.

But is there any simpler method?

-----advanced version

what if I want X and Y and Z all be AR(1) process? I am running EnKF by perturbing atmospheric forcing. Papers just said they made perturbations like this, but I don't know how do it. The data is hourly precipitation, however they said "The time series correlation (temporal correlation) is applied to a first-order autoregressive model." Table 2 shows that they apply a temporal correlation of 24 hours...(https://doi.org/10.1175/JHM-D-15-0037.1)

Thanks!

Answer 1

Strategy:

1. Calculate the normal mean and variance for the lognormal variables, then simulate the normal variables and calculate $e^\cdot$.
2. An AR-process with diagonal $\phi$ and noise cov $C_a$ will have variance $C_a / (1-e)$ where $e$ is given below. Therefore, simulate the AR-process with target cov $C(1-e)$ so that the realized AR-cov becomes $C$.

import numpy as np
from numpy import diag, log, exp, ones
G = 10**6
ms  = [1,1,0]
s   = [0.5,0.3,20]
mx  = ms[0]
my = ms[1]
rho = np.array([[1, -0.8, 0.5], [-0.8, 1, -0.5], [0.5, -0.5, 1] ])
C   = diag(s) @ rho @ diag(s)   # Target cov
# AR model: a' = phi a + U[t]
phi = np.array([0.2, 0.30, 0.2]) # Set this to zero for non-ar case
# (a',a)_xy = (phi a_x + X, phi a_y + Y) = [..] = C[X,Y]/(1-phi_x*phi_y)
e = diag(phi) @ ones([len(phi)]*2) @ diag(phi)
Ca = C * (1 - e)      # Target cov for AR
s2w = log(Ca[0,0] / mx**2 +1) # V[W] = V[log(X)]
s2v = log(Ca[1,1] / my**2 +1) # V[V] = V[log(Y)]
m = ms.copy()                 # new mean
m[0] = log(mx) - s2w/2        # E[W] = E[log(X)]
m[1] = log(my) - s2v/2        # E[V] = E[log(Y)]
CC = Ca.copy()           # new cov
CC[:,0] = CC[:,0] / mx   # C[X,Z] = C[W,Z] * mx so div by mx gives new cov
CC[0,:] = CC[0,:] / mx
CC[0,0] = log(CC[0,0] + 1)
CC[:,1] = CC[:,1] / my   # C[Y,W] = C[V,W] * my so div by my gives new cov
CC[1,:] = CC[1,:] / my
CC[1,1] = log(CC[1,1] + 1)
CC[0,1] = CC[1,0] = Ca[1,0] / mx / my    # C[X,Y] = C[W,Y] * mx = C[W,V] * mx * my
U = np.random.randn(G,len(m)) @ np.linalg.cholesky(CC).T + m
U[:,0] = exp(U[:,0])  # if phi = 0, U has mean ms and cov C
U[:,1] = exp(U[:,1])  # if phi = 0, U has mean ms and cov C
Part 2: Sim AR
a = np.zeros((G,len(m))) # sim AR:
for t in range(1, G): a[t] = phi * a[t-1] + (U[t] - ms)
a = a + ms # add mean
print('    Simulated Cov / target Cov')
print(np.cov(a.T) / C)
print('    Simulated mean - target mean')
print(a.mean(0) - ms)

The parameters impose restrictions on $\phi$. They must take values so that $C(1 - e(\phi))$ is positive definite.

Edit: Y is now also log normal per the posters request.

If you, for instance, want a temporal correlation such that $Z$ 24 hours ago is correlated to $X,Y$, then you could simply take the simulation above and lag $Z$:

L = 24
ag = a.copy()
ag[:-L, 2] = ag[L:, 2]   # Send simulated z back in time
print(np.corrcoef(ag[L:,0], ag[:-L,2])[1,0]) # corr(ax[t], az[t-24])

This has zero correlation on all horizons except L=24. For better control, and perhaps a more natural decay of correlation, you can specify $\phi$ as a matrix. This makes things slightly more complicated and requires you to calculate the $\phi$-matrix needed for the desired correlation.

The value for $\phi$ depends on what you are trying to simulate. I would make a rough estimate (either by looking at data or numerically estimating) set $\phi$ to similar values. If these values are creating non-PSD Covs, you would have to change your assumptions about the covariance.

Answer 2

To generate random numbers from correlated distributions where two are lognormal and one is normal, and then extend it to an AR(1) process, you can follow these steps:

Step 1: Define Parameters and Correlation Matrix

First, set up the means, standard deviations, and the correlation matrix.

import numpy as np
from scipy.stats import multivariate_normal
Define the means and standard deviations
mean = [1, 1, 0]
std_dev = [0.5, 0.3, 20]
Define the correlation matrix
corr_matrix = np.array([[1, -0.8, 0.5],
                        [-0.8, 1, -0.5],
                        [0.5, -0.5, 1]])

Calculate the covariance matrix from the correlation matrix

cov_matrix = np.diag(std_dev).dot(corr_matrix).dot(np.diag(std_dev))

Step 2: Generate Multivariate Normal Distribution

Generate a multivariate normal distribution using these parameters. This will be transformed for the lognormal distributions.

# Generate samples from multivariate normal distribution
mvn_samples = multivariate_normal(mean, cov_matrix).rvs(size=1000)

Step 3: Transform to Lognormal Distributions for X and Y

Transform the first two variables into lognormal distributions.

# Convert to lognormal for X and Y
mvn_samples[:, 0] = np.exp(mvn_samples[:, 0])
mvn_samples[:, 1] = np.exp(mvn_samples[:, 1])
X, Y, Z now contain the desired distributions
X, Y, Z = mvn_samples.T

Step 4: Generate AR(1) Processes for X, Y, and Z

Define a function to generate an AR(1) process and apply it to X, Y, and Z.

def generate_ar1_process(alpha, noise_series):
    series = [noise_series[0]]
    for i in range(1, len(noise_series)):
        series.append(alpha * series[i-1] + noise_series[i])
    return series
alpha = 0.8  # Autoregression coefficient
X_ar1 = generate_ar1_process(alpha, X)
Y_ar1 = generate_ar1_process(alpha, Y)
Z_ar1 = generate_ar1_process(alpha, Z)

Visualize the distributions and the AR(1) processes.

import matplotlib.pyplot as plt
Plotting the transformed distributions
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.hist(X, bins=30, color='blue', alpha=0.7)
plt.title('Transformed X (Lognormal)')
plt.subplot(1, 3, 2)
plt.hist(Y, bins=30, color='green', alpha=0.7)
plt.title('Transformed Y (Lognormal)')
plt.subplot(1, 3, 3)
plt.hist(Z, bins=30, color='red', alpha=0.7)
plt.title('Unchanged Z (Normal)')
plt.tight_layout()
plt.show()

# Visualization of the AR(1) Processes
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.plot(X_ar1, color='blue')
plt.title('AR(1) Process for X')
plt.subplot(1, 3, 2)
plt.plot(Y_ar1, color='green')
plt.title('AR(1) Process for Y')
plt.subplot(1, 3, 3)
plt.plot(Z_ar1, color='red')
plt.title('AR(1) Process for Z')
plt.tight_layout()
plt.show()

Note

I'm not 100% sure I got all your requirements right. An optional step you could add at the start of the code is to adjust the mean and covariance matrix calculations. You could do it like so:

def lognormal_params(mean, std_dev):
    var = std_dev ** 2
    mean_log = np.log(mean ** 2 / np.sqrt(var + mean ** 2))
    std_dev_log = np.sqrt(np.log(var / mean ** 2 + 1))
    return mean_log, std_dev_log
Mean and std dev for X and Y (lognormal) and Z (normal)
mean_X, std_X = 1, 0.5
mean_Y, std_Y = 1, 0.3
mean_Z, std_Z = 0, 20
Adjusted mean and std dev for underlying normal distributions of X and Y
mean_X_log, std_X_log = lognormal_params(mean_X, std_X)
mean_Y_log, std_Y_log = lognormal_params(mean_Y, std_Y)
Updated means and standard deviations
mean = [mean_X_log, mean_Y_log, mean_Z]
std_dev = [std_X_log, std_Y_log, std_Z]
...

Hopefully it gives you some material to solve your problem.

Answer 3

Based on @Hunaphu's answer, I made a light change and here is a modified version:

import numpy as np
from numpy import diag, log, exp, ones
G = 10**6
ms  = [1,1,0]
s   = [0.5,0.3,20]
mx  = ms[0]
my = ms[1]
rho = np.array([[1, -0.8, 0.5], [-0.8, 1, -0.5], [0.5, -0.5, 1] ])
C   = diag(s) @ rho @ diag(s)   # Target cov
# AR model: a' = phi a + U[t]
phi = np.array([0.9, 0.9, 0.9]) # Set this to zero for non-ar case
# (a',a)_xy = (phi a_x + X, phi a_y + Y) = [..] = C[X,Y]/(1-phi_x*phi_y)
e = diag(phi) @ ones([len(phi)]*2) @ diag(phi)
# Ca = C * (1 - e)      # Target cov for AR
Ca = C * (1)      # Target cov for AR
s2w = log(Ca[0,0] / mx**2 +1) # V[W] = V[log(X)]
s2v = log(Ca[1,1] / my**2 +1) # V[V] = V[log(Y)]
m = ms.copy()                 # new mean
m[0] = log(mx) - s2w/2        # E[W] = E[log(X)]
m[1] = log(my) - s2v/2        # E[V] = E[log(Y)]
CC = Ca.copy()           # new cov
CC[:,0] = CC[:,0] / mx   # C[X,Z] = C[W,Z] * mx so div by mx gives new cov
CC[0,:] = CC[0,:] / mx
CC[0,0] = log(CC[0,0] + 1)
CC[:,1] = CC[:,1] / my   # C[Y,W] = C[V,W] * my so div by my gives new cov
CC[1,:] = CC[1,:] / my
CC[1,1] = log(CC[1,1] + 1)
CC[0,1] = CC[1,0] = Ca[1,0] / mx / my    # C[X,Y] = C[W,Y] * mx = C[W,V] * mx * my
CCa = CC * (1 - e)      # Target cov for AR
U = np.random.randn(G,len(m)) @ np.linalg.cholesky(CCa).T + m* (1 - phi)
Part 2: Sim AR
a = np.zeros((G,len(m))) # sim AR:
for t in range(1, G): a[t] = phi * a[t-1] + (U[t])
a_tmp = a
a_tmp[:,0] = np.exp(a[:,0])
a_tmp[:,1] = np.exp(a[:,1])
print(np.min(a_tmp, axis=0))
print('    Simulated Cov / target Cov')
print(np.cov(a_tmp.T) / C)
print('    Simulated mean - target mean')
print(a_tmp.mean(0) - ms)

So basically define the AR1 process of lognormal distribution as:$\ln Y_{t} = \phi * \ln Y_{t-1}+N(\mu, \sigma^2)$ where the $Y_{t}$ follows a lognormal distribution and the $N(\mu, \sigma^2)$ is a normal distribution. So our strategy is to sample from a multivariate normal distribution, then convert it to an AR1 process, finally use np.exp to make it lognormally distributed. The core part is to generate random noise of $N(\mu, \sigma^2)$ to satisfy the certain values of the mean and std of the $Y_{t}$. There are two transformations:

The first one is to convert $\mu_{normal}$ and $std_{normal}$ to the $\mu_{lognormal}$ and $std_{lognormal}$. Then using $\mu_{lognormal}$ and $std_{lognormal}$ we can get the mean and std for $N(\mu, \sigma^2)$.