Create random time series with shifts in R I'm trying to create some random time series with shifts along the series.
Some with normal distribution and the others with non-normal distribution (Log normal for example).
For both I have to decide:


*

*the numbers of the shifts

*the positions of each shifts

*the magnitudes of each shifts

*the length between a shift and another one


The data is 130 years length (monthly) with all positive numbers (because I want to simulate natural abundance time series).
I hope with a picture you can understand better.

I've foun in other questions some methods for example arima.sim. But I do not think that I can put in my shifts.
 A: You will have to generate the series by yourself. Fortunately it is not that hard. Simple ARMA(p,q) model:
$$Y_t=\alpha_1 Y_{t-1}+...+\alpha_pY_{t-p}+Z_t+\theta_1Z_{t_1}+...+\theta_qZ_{t-q}$$
can be generated with simple loop:
Z<-rnorm(1000)
Y<-rep(NA,1000)
a <- runif(p)
b <- runif(q)

for(i in (max(p,q)+1):1000) {
     Y[i] <- a*Y[i-(1:p)]+c(1,b)*Z[i-(0:q)]
}

It is advisable to use a burn-in period, i.e. generate additional $n_b$ values and discard the first $n_b$, so that the series should not depend on the initial values.
Now your shifts can be a simple vector $S$ the same as the length of $Y$ with $i$-th element containing the magnitude of the shift, if there is one, and zero otherwise. Then the loop can be modified as following:
for(i in (max(p,q)+1):1000) {
     Y[i] <- S[i]+ a*Y[i-(1:p)]+c(1,b)*Z[i-(0:q)]
}

$S$ can be generated as follows:
## Number of shifts:
ns <- rpois(1, 10)
## The length of shifts
ls <- rpois(ns, 10)
## The positions of the shifts
ps <- cumsum(ls)
## The magnitudes of the shifts
ms <- rchisq(ns, 10)
##Shift vector
S <- rep(0, length(Y))
S[ps] <- ms

Naturally you need to choose your own distributions, but the idea is the same. 
In R loops can be slow, but it is easy to speed up them in this case with function filter. Look at the examples in the help page or examine the code of arima.sim for ideas. 
Note that for ARMA(p,q) processes the coefficients must satisfy certain constraints. I suggest then to start with ARMA(1,1) model with $|\alpha_1|<1$ and $|\theta_1|<1$.
