;+
; NAME:
;	MONTE_CARLO_AR1
;
; PURPOSE:
;	This function returns a Monte Carlo generated time series which 
;	preserves the distributional and AR(1) properties of the original 
;	series.
;
; CATEGORY:
;	Time Series Analysis
;
; CALLING SEQUENCE:
;	Result = MONTE_CARLO_AR1( X )
;
; INPUTS:
;	X:  An array of any numerical type.
;
; KEYWORD PARAMETERS:
;	REPLACE:  If set, the function shuffles with replacement.  The default
;		is without replacement.
;	SEED:  A seed for the random number generator to use.
;
; OUTPUTS:
;	Result:  Returns a Monte Carlo generated vector inspired by X.
;
; USES:
;	SHUFFLE.pro
;
; PROCEDURE:
;	This function generates an AR(1) vector with the random component 
;	generated by shuffling the elements of the input vector.
;
; EXAMPLE:
;	Define an AR(1) vector with coefficient 0.9, length 1000, and unit 
;	standard deviation noise.
;         x = fltarr( 1000 )
;         for i = 1, 999 do x[i] = 0.9 * x[i-1] + randomn( seed, 1 )
;	Generate an time series that maintains the distribution and AR(1) 
;	properties of x.
;	  result = monte_carlo_ar1( x, replace=1 )
;
; MODIFICATION HISTORY:
; 	Written by:	Daithi A. Stone (stoned@atm.ox.ac.uk), 2003-11-19.
;	Modified:	DAS, 2005-10-23 (added SEED keyword)
;-

;***********************************************************************

FUNCTION MONTE_CARLO_AR1, $
	X, $
	REPLACE=replaceopt, $
	SEED=seed

;***********************************************************************
; Constants and Variables

; Length of input vector
nx = n_elements( x )

; Estimate the AR(1) regression coefficient
alpha1 = correlate( x[0:nx-2], x[1:nx-1] )

; The mean of the input vector
xmean = mean( x )

;***********************************************************************
; Generate New Series

; Initialise new vector
y = fltarr( nx )

; Shuffle the entries in the input vector
xshuf = shuffle( x-xmean, replace=replaceopt, seed=seed )

; Pick first entry in new vector.
; The factor ensures the variance and its partitioning are identical for the 
; input and new vectors
y[0] = xshuf[0] * sqrt( 1. - alpha1^2 )

; Pick the rest of the entry in the new vector.
; The factor ensures the variance and its partitioning are identical for the 
; input and new vectors
for i = 1, nx - 1 do y[i] = alpha1 * y[i-1] + xshuf[i] * sqrt( 1. - alpha1^2 )

; Restore the mean to the new vector
y = y + xmean

;***********************************************************************
; The End

return, y
END