We're looking to write a GAMS parameter that generates a series of random normal numbers with a user-chosen mean and user-chosen variance. But here's the twist: We need this generated series of numbers to have a user-chosenauto-correlationfrom one observation to the next. So if one number in the sequence exceeds the mean, we need the next one to be above it, too.

We're actually thinking of streamflows from one month or year to the next. We know the flows' means and variances but also know their autocorrelation. Do you know of any kind of GAMS function that could generate this kind of parameter. Possibly you could give us tips on how to write this function in GAMS?

I suspect you want to simulate an AR(1) process here. I.e.

X(t) = φ

_{0}+ φ_{1}X(t-1) + ε(t)

where ε(t) is normally distributed with mean 0 and variance σ_{e}^{2}.

Now we have the following:

- Mean of x(t): μ = φ
_{0}/(1−φ_{1}) - Variance of x(t): σ
^{2}=σ_{e}^{2}/(1−φ_{1}^{2}) - Autocorrelation coefficient: ρ=φ
_{1}

So given (μ,σ^{2},ρ) we can calculate

- φ
_{1}=ρ - φ
_{0}= μ(1−φ_{1}) - σ
_{e}^{2}=σ^{2}(1−φ_{1}^{2})

Now we can simulate X(t). In the hope that I did all the math correctly, we can now do in GAMS:

$ontext

simulate ar(1) process

x(t) = phi0 + phi1*x(t-1) + N(0,esigma2)

mean of x(t): mu = phi0/(1-phi1)

variance of x(t): var = esigma2/(1-phi1^2)

autocorrelation: rho = phi1

==>

given mu,var,rho we have

phi1 = rho

phi0 = mu*(1-phi1)

esigma2 = var*(1-phi1^2)

$offtext

set

i0 'warmup-period' /i0-1*i0-100/

i 'number of observations to simulate' /i1*i50/

;

** these are the inputs* mu = mean of x(t)* s2 = variance of x(t)* rho = autocorrelation coefficientscalar

mu 'mean' /0/

s2 'variance' /2/

rho 'autocorrelation' /0.8/

;

scalars

phi0

phi1

esigma

;

phi1 = rho;

phi0 = mu*(1-phi1);

esigma = sqrt(s2*(1-sqr(phi1)));

** warm up*scalarx0 /0/;loop(i0,

x0 = phi1*x0 + normal(phi0,esigma);

);

** simulate AR(1) process*parameterx(i);loop(i,

x0 = phi1*x0 + normal(phi0,esigma);

x(i) = x0;

);

displayx;

You may want to throw x(i) into a statistics package and see if it really behaves as you want.

Note: Instead of the warm-up loop, we could also start with a first observation from N(μ,σ^{2}).

(Compare normal with AR1/rho=0.8: the AR(1) series stays positive or negative longer)