How to create noise with time dependent correlation in Julia?

I want to create real noise $\xi(t)$ which fulfills the property $E[xi_t\xi_s] = f(t, s)$ in Julia. Is there a library or function with which I could implement this?
I tried the impementation with $f(t, s) = e^{-(t-s)}$
My code is not working (the noise is decreasing exponentially).

using LinearAlgebra

function generate_gaussian_noise(timesteps)
    n = length(timesteps)
    noise = zeros(n)

    for i in 1:n
        t = timesteps[i]
        noise[i] = real(randn())
    end

    return noise
end

function construct_correlated_noise(timesteps)
    n = length(timesteps)
    noise = generate_gaussian_noise(timesteps)
    correlation_matrix = zeros(n, n)

    for i in 1:n
        for j in 1:n
            correlation_matrix[i, j] = exp(-(timesteps[i] - timesteps[j]))
        end
    end

    # Perform eigenvalue decomposition
    eigenvalues, eigenvectors = eigen(correlation_matrix)
    D = Diagonal(sqrt.(eigenvalues))
    Q = eigenvectors
    L = Q * D

    correlated_noise = L * noise

    return correlated_noise
end

# Example usage
timesteps = collect(1:401)
correlated_noise = construct_correlated_noise(timesteps)

If I understand you correctly you want something like this?

using Distributions
f(t, s) = exp(-abs(t-s))
timesteps = 1:401;
covariance_matrix = f.(timesteps, timesteps');
mean_vector = zeros(length(timesteps));
d = MvNormal(mean_vector, covariance_matrix);
x = rand(d) # this is a single sample

Note that f(t,s) must have a property that f(t,s)=f(s,t) which your function did not have, that is why I have added abs. I also modelled covariance because what you defined is not correlation but covariance assuming that the expected value of the process is constant and equal to 0.

It looks like you're just asking how to simulate a zero-mean Gaussian vector with arbitrary correlation structure. You can just use the Distributions.jl package to do this. For example:

using Distributions
numobs = 100;
muvec = zeros(Float64, numobs);
covmat = Matrix{Float64}(I, numobs, numobs);
x = rand(MvNormal(muvec, covmat);

You will, of course, need to populate the covmat line with whatever covariance matrix structure you are after. I've just used code for a diagonal covariance matrix here, but you can do whatever you want.

If you have a more specific correlation structure in mind, e.g. maybe an autoregressive structure, then there will be more efficient ways to go about this. But for an arbitrary structure, this is probably the best and easiest approach.