Is this AR(1)?

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-chosen auto-correlation from 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) = φ₀ + φ₁ X(t-1) + ε(t)

where ε(t) is normally distributed with mean 0 and variance σ_e².

Now we have the following:

• Mean of x(t): μ = φ₀/(1−φ₁)
• Variance of x(t): σ²=σ_e²/(1−φ₁²)
• Autocorrelation coefficient: ρ=φ₁

So given (μ,σ²,ρ) we can calculate

• φ₁=ρ
• φ₀= μ(1−φ₁)
• σ_e²=σ²(1−φ₁²)

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 coefficient
scalar
   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
*
scalar x0 /0/;
loop(i0,
   x0 = phi1*x0 + normal(phi0,esigma);
);

*
* simulate AR(1) process
*
parameter x(i);
loop(i,
   x0 = phi1*x0 + normal(phi0,esigma);
   x(i) = x0;
);

display x;

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(μ,σ²).

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