Time Series Seasonality (cycle) not constant I'm having problems creating a reliable model to forecast this time series (quarterly data):
3235 3050 3045 2507 2346 2360 2396 2577 2518 2518 2683 2714 2711 2711 2731 2551 2474 2516 2526 3264 2569 2480 2885 3028 2913 3092 3091 3202 3143 3256 3106 2917 3184 3381 3343 3608 4152 3948 3839 3474 

I'm using the library(forecast) from R. I specified several models like this:
model1 <- Arima(mydata, order=c(0,1,0), seasonal=list(order=c(0,1,1),period=4))
model2 <- Arima(mydata, order=c(0,0,0), seasonal=list(order=c(0,1,1),period=4),include.constant = TRUE)

Variations of the last one with mean and drift and many more.
I get poor results when I test the models with 80% of the data over the remaining 20%. 
The issue is the seasonality (cycle) because my data begins with one peak, then in the 20th period has another peak and then in periods 37th, 38th and 39th peaks again. So the seasonality is not constant every 20 periods and based on latest data points I can have more than peak next to the other. I changed the parameter "period=4" to "period=17 or 18 or 19", but results still far from optimal.   
I tried with external variables without success. 
Any ideas or advice?
Should I be using something different from ARIMA? Maybe via intervention variables? (I don't know how to do it in R). Any smoothing technique? 
Can anybody provide some code? Or point out to methods or libraries for this problem?
 A: You have 40 data points, and it is quarterly data, your seasonality should be yearly.
According to this plot, I don't see any evidence for seasonality. There is possible an additive trend component at best.

A: Your 40 values do exhibit a quarterly effect requiring seasonal differences. You need to incorporate interventions. 
The ACF of the original series is dampened by the 5 exceptional values/pulses/intereventions which tend to downwards bias the true acf due to the inflated variance.  . 
The AIC should be calculated from residuals using models that control for intervention administration, otherwise the intervention effects are taken to be Gaussian noise, underestimating the actual model's autoregressive effect and thus miscalculates the model parameters which leads directly to an incorrect error sum of squares and ultimately an incorrect AIC and ultimately bad model identification.
Additionally there are 5 identifiable pulses . and  and  leading to the following acf of the residuals  suggesting sufficiency. The model is (1,0,0)(0,1,0) 4 .
