How to generate a standard normal distributed time-series with a given ACF

Question

I want to generate a standard normal distributed time-series. In addition the ACF of my timeseries should match a desired ACF. I have given lags 1 to 30 with the corresponding ACF-values. For further processing, the CDF of the timeseries should be the standard normal CDF.

When i try filtering a standard normal distributed timeseries with my ACF, the result does not have the desired attributes. Neither is the ACF like the desired ACF, nor is the CDF unchanged.

I tried the approach from jblood94 from this and this thread. But this approach only works for very specific ACF-values.

When i insert my desired ACF, this approach does not work.

Here is my desired ACF:

The model for my ACF:

lags   =  0:30
a      =  0.9999
b      = -0.07197 
acf    = a*exp(b.*lags)

I assume that the problem is the high correlation at the beginning.

The timeseries needs to be standard normal distributed, for the next steps in my program.

The following Matlab-Code includes the approach. But as you can see, it will only work, if we set a minimum of 45 lags. Why is this so? In addition, i have to adjust the initial STD of my samples. How can i calculate the initial STD, to be after the filtering 1?

%% init
clear;
close all;
clc;
set(groot, 'DefaultAxesLineWidth', 2);
set(groot, 'DefaultLineLineWidth', 2);
set(groot, 'DefaultAxesFontSize', 14);
%% define desired auto-correlation function
lagsDesired = 0:45; % if lag is <45 the approach fails
a =  0.999;
b = -0.07197;
acfDesired = a * exp(b .* lagsDesired);
%% generate normal distributed samples 
numSamples = 1e6;
% sigma increases during filtering (How can i calculate the factor?)
sigmaFactor = 0.35;
x = sigmaFactor * randn(numSamples, 1);
%% predefine indices
m = length(acfDesired);
i1 = sequence((m-1):-1:1);
i2 = sequence((m-1):-1:1, 'from', 2:m);
i3 = cumsum((m-1):-1:1);
%% initial values for adjustment
w = acfDesired;
temp = cumsum(w(i1) .* w(i2));
a = [1, diff([0  temp(i3)])/sum(w.^2)];
%% iterative adjustment of the filter weights
while (max(abs(acfDesired - a)) >= 2e-3) %sqrt(eps))
    max(abs(acfDesired - a))
    temp = cumsum(w(i1) .* w(i2));
    a = [1, diff([0  temp(i3)])/sum(w.^2)];
    w = w .* (acfDesired./a);
end
%% filter normal samples with weights
y = filter(w, 1, x);
%% transform normal to uniform samples via standard normal cdf
uniRandCorr = normcdf(y);
%% plot
figure; 
[Fy, xy] = ecdf(y);
[Fn, xn] = ecdf(randn(numSamples, 1));
plot(xy, Fy);
hold on; 
plot(xn, Fn);
l = legend({'ECDF of samples' , 'Desired ECDF'}, 'Location', 'best');
title(l, 'Legend');
title('CDF');
xlabel('Value');
ylabel('CDF');
grid on;
figure;
[acf, lags] = autocorr(uniRandCorr, 'NumLags', 60);
plot(lags, acf);
hold on; 
plot(lagsDesired, acfDesired);
l = legend({'ACF of samples' , 'Desired ACF'}, 'Location', 'best');
title(l, 'Legend');
title('ACF');
xlabel('Lag');
ylabel('ACF');
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [seq] = sequence(nvec, varargin)
% [seq] = sequence(nvec, varargin)
%
% DESCRIPTION:
%   R-Functions translated to MATLAB
%   Creates a vector of sequences
%   R: A Language and Environment for Statistical Computing Reference Index
%   Version 4.2.3 Page 505 ff
%
% EXAMPLE:
%   -
%
% INPUT:
%   nvec                    :   Input vector
%   from                    :   Starting value(s)
%   by                      :   Increment of the sequence(s)
%
% OUTPUT:
%   seq                     :   Generated sequence
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parse inputs
iParser = inputParser;
addRequired(iParser, 'nvec', ...
    @(x) isnumeric(x) && all(x>=0));
addOptional(iParser, 'from', ones(length(nvec), 1), ...
    @(x) isnumeric(x) && (isscalar(x) || length(x)==length(nvec)));
addOptional(iParser, 'by', ones(length(nvec), 1), ...
    @(x) isnumeric(x) && (isscalar(x) || length(x)==length(nvec)));
parse(iParser, nvec, varargin{:});
if isscalar(iParser.Results.from)
    from = ones(length(nvec))*iParser.Results.from;
else
    from = iParser.Results.from;
end
if isscalar(iParser.Results.by)
    by = ones(length(nvec))*iParser.Results.by;
else
    by = iParser.Results.by;
end
seq = [];
for i=1:length(nvec)
    seq = [seq, linspace(from(i), from(i) + (by(i)*(nvec(i)-1)), nvec(i))];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This gives the following output:

As you can see, the result is correct. BUT, it works only with min. 45 lags AND if i change the model of the acf I have to adapt this parameter. In addition i have no idea, how i can calculate the sigma of the gauss samples, so that the sigma AFTER the filtering is 1 (standard). I know, that the standarddeviation changes during filtering, but I have not found a formula for this yet.

Answer 1

Theoretical Backgorund

For the sake of simplicity let's assume $\sigma^2=\gamma(0)=1$, where $\gamma(h)$ represents the autocovariance at lag $h$. Let your time series be $\{X_i\}_{i=1}^n$. Consider the covariance matrix of your data

$$\begin{align} \Sigma_X=\begin{bmatrix}1 & \gamma(1) & \gamma(2) & \cdots & \gamma(n-1)\\\gamma(1) & 1 & \gamma(1) & \cdots & \gamma(n-2)\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \gamma(n-1) & \gamma(n-2) & \gamma(n-3) & \cdots & 1\end{bmatrix}\end{align}.$$

For a stationary process (I am assuming you would like your time series to be stationary) this covariance matrix will be positive semi-definite (see https://stat.ethz.ch/Teaching/buhlmann/Zeitreihen/acf.pdf, and Show that the autocovariance function of stationary process {${X_t}$} is positive definite). Now since $\Sigma_X$ is positive semi-definite and symmetric, it will have a symmetric square root such that $$\Sigma_X=\Sigma_X^{1/2}\Sigma_X^{1/2}.$$

Finally we can use the fact that for a vector of random variables $Z$, $\text{Cov}(AZ)=A\Sigma_ZA^T$. So here if $Z$ is a vector of normal random varaibles with $\Sigma_Z=I$ then $$\text{Cov}(\Sigma_X^{1/2}Z)=\Sigma_X.$$

Application in R

So in your case what you would do is

1. Using the ACF you have calculated, create a $\Sigma_X$ matrix in R.
2. Calculate $\Sigma_X^{1/2}$ (you could use the sqrtm() function in R).
3. Generate $Z$ using rnorm() in R, and then multiply $Z$ by $\Sigma_X^{1/2}$ to get your time series with the desired covariance, and hence autocovariance, structure.

Answer 2

acf    = a*exp(b.*lags)

You can get this type of autocovariance function with an AR1 process $$X_{t} = \phi X_{t-1} + w_{t} $$

where the $w_t$ are independent and identical distributed as $N(0,\sigma^2)$

this gives an acf of the form:

$$\text{ACF}(\text{lag}) = \frac{\sigma^2}{1-\phi^2} \phi^{\text{lag}} $$

(see the question Finding the ACF of AR(1) process)

You can find $\phi$ and $\sigma$ from your $a$ and $b$ if you compare with your function and set $\text{lag} = 0$ and $\text{lag} = 1$

$$\frac{\sigma^2}{1-\phi^2} = a$$ $$\frac{\sigma^2}{1-\phi^2} \phi = a \exp(-b)$$

from which it follows that $\phi = -\exp(b)$ and $\sigma^2 = a[1-\exp(-2b)]$.