ar1etp: AR(1) + ETP simulation routine

Description

Simulate a combined AR(1) + ETP signal, plot spectrum and confidence levels

Usage

ar1etp(etpdat=NULL,nsim=100,rho=0.9,wtAR=1,sig=90,tbw=2,padfac=5,ftest=F,fmax=0.1,
       speed=0.5,pl=2,graphfile=0)

Arguments

etpdat: Eccentricity, tilt, precession astronmical series. First column = time, second column = ETP. If not entered will use default series from Laskar et al. (2004), spanning 0-1000 kyr.
nsim: Number of simulations.
rho: AR(1) coefficient for noise modeling.
wtAR: Multiplicative factor for AR1 noise (1= eqivalent to ETP variance). If < 0, etp signal will be excluded from the simulations (noise only)
sig: Demarcate what confidence level (percent) on plots?
tbw: MTM time-bandwidth product.
padfac: Pad with zeros to (padfac*npts) points, where npts is the number of data points.
ftest: Include MTM harmonic f-test results? (T or F)
fmax: Maximum frequency for plotting.
speed: Set the amount of time to pause before plotting new graph, in seconds.
pl: Plot (1) log frequency-log power or (2) linear frequency-linear power?
graphfile: Output a pdf or jpg image of each plot? 0 = no, 1 = pdf, 2 = jpeg. If yes, there will be no output to screen. Individual graphic files will be produced for each simluation, for assembling into a movie.

Details

Note: Setting wtAR=1 will provide equal variance contributions from the etp model and the ar1 model. More generally, set wtAR to the square root of the desired variance contribution (wtAR=0.5 will generate an AR1 model with variance that is 25% of the etp model). If you would like to exclusively evaluate the noise (no etp), set wtAR < 0.

Note: You may use the function etp to generate eccentricity-tilt-precession models.

References

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A.C.M., Levrard, B., 2004, A long term numerical solution for the insolation quantities of the Earth: Astron. Astrophys., Volume 428, 261-285.

Examples

Run this code

 if (FALSE) {

# run simulations using the default settings
ar1etp()

# compare with a second model:
# generate etp model spanning 0-2000 ka, with sampling interval of 5 ka.
ex1=etp(tmin=0,tmax=2000,dt=5)
# run simulations, with rho=-.7, and scaling noise to have 50
ar1etp(etpdat=ex1,rho=0.7,wtAR=sqrt(0.5))

 }