Solved – Analysing some series in R – problems with TS/decompose functions

decompositionforecastingfrequencyrtime series

I am trying to analyze and forecast the European Emission Allowance (EUA) prices by using the data from January 2008 to March, 2016 and also using other series that might be related to that price, as Brent or IBEX35 prices (I will use correlations or divergences between those series).

The problem is that, in my case, I do not know how to modify the frequency correctly when using the ‘ts’ function on my series because I have taken the data from investing.com and sendeco2.com and the time frequency they use is just everyday the stock market opens, so ‘365’ should not be the correct frequency value… I have checked how many values I have for each year and the results are:

-2008: 252 values
-2009: 250 values
-2010: 252 values
-2011: 249 values
-2012: 254 values
-2013: 254 values
-2014: 255 values
-2015: 254 values
-2016 (March): 105 values

It means that my series are irregular. What i do is to store each series (EUA, Brent…etc) to a vector (total length of 2125 for each one) and I work with those variables:

BDD<-read.csv('BDD.csv',colClasses=c("Date","numeric","numeric","numeric","numeric"), header=TRUE, sep=";")
Date <- BDD[1:nrow(BDD),1];
EUA <- BDD[1:nrow(BDD),2]; 
Brent <- BDD[1:nrow(BDD),3];
IBEX35 <- BDD[1:nrow(BDD),4];
PME <- BDD[1:nrow(BDD),5];

I have also tested with freq=1 (but it is obviously a mistake because they are not yearly values) and i got this error:

> EUA_ts<-ts(EUA, frequency=1, start=c(2008,1), end=c(2016,3))
> EUA_decomp<-decompose(EUA_ts, type=c("additive"))
Error in decompose(EUA_ts, type = c("additive")) : 
  time series has no or less than 2 periods

Is there another possibility to decompose a series? Am I doing something wrong?

I have uploaded my database in Drive to bring out more detail, if necessary:
https://drive.google.com/open?id=0B7nP03_LfDvQZS1WNUdqbWM1X1U

Best Answer

There is an implied statement in your question that I want to clarify before giving an answer. You stated :

I have also tested with freq=1 (but it is obviously a mistake because they are not yearly values)

This implies that either data frequency is only governed by frequency per year or that all and or only yearly data have a frequency of 1. Both of these statements are false. Solar cycles have a frequency of about 11 years, so yearly solar cycle data would have a frequency of about 11. On the other hand, something like server copy errors per copy attempt may be sampled more than once per second, but still have a frequency of 1 since copy errors are a mostly random process. When determining the frequency of the series, you need to think about what is causing the seasonality and not just rely on the sample rate (though that can be a good starting point). The answer below assumes you have correctly identifies the seasonal period as 1 year.

As to your original question, there are a few ways you can deal with your data. The simplest would be to set your seasonality to the mean value of your number of samples per year. Time series do not need to have an integer seasonality, so you can make a time series in R as below:

ts(rnorm(100), frequency = 14.73)

Another option would be to add back in the missing days to the time series and use a frequency of 365 (or even better 365.24). If you have the time points for each observation, you can use the zoo package in R to make an irregular time series and then fill in the missing values.

You can make the series using:

x.Date <- as.Date("2003-02-01") + sample(1000,900) - 1
x <- zoo(rnorm(900), x.Date)
y <- ts(as.ts(x), frequency = 365.24)

The missing values can be filled using many methods, but one to consider is zoo::na.approx. From there the series can be decomposed as normal.

decompose(na.approx(y))

A few final notes and options: 1) decompose is a useful method, but you may also want to consider the stl decomposition for your data as even the decompose docs say that "stl provides a much more sophisticated decomposition." 2) You will always get an error if you try to use decompose or stl on a series with a frequency of 1. Both functions seek to separate the seasonal and trend components of the data, so if there is no seasonal component (ie frequency = 1), there is a problem. If instead you just want to separate trend from noise, you might consider using a moving average.

Related Solutions

Solved – How to decompose a time series with multiple seasonal components

R's forecast package bats() and tbats() functions can fit BATS and TBATS models to the data. The functions return lists with a class attribute either "bats" or "tbats". One of the elements on this list is a time series of state vectors, $x(t)$, for each time, $t$.

See http://robjhyndman.com/papers/complex-seasonality/ for the formula's and Hyndman et al (2008) for a better description of ETS models. BATS and TBATS are an extension of ETS.

For example:

fit <- bats(myTimeseries)
fit$x

In this case, each row of x will be on fourier-like harmonic.

There are also plot.tbats() and plot.bats() functions to automatically decompose and view the components.

R – Solving Problems with Time Series Prediction

You are attempting to forecast a compositional time series. That is, you have three components that are all constrained to lie between 0 and 1 and to add up to 1.

You can address this issue using standard exponential smoothing, by using an appropriate generalized logistic transformation. There was a presentation on this by Koehler, Snyder, Ord & Beaumont at the 2010 International Symposium on Forecasting, which turned into a paper (Snyder et al., 2017, International Journal of Forecasting).

Let's walk though this with your data. Read the data into a matrix obs of time series:

obs <- structure(c(0.03333333, 0.03810624, 0, 0.03776683, 0.06606607, 
0.03900325, 0.03125, 0, 0.04927885, 0.0610687, 0.03846154, 0, 
0.06028636, 0.09646302, 0.04444444, 0.01111111, 0.02309469, 0.03846154, 
0.03119869, 0.01201201, 0.02058505, 0.015625, 0, 0.01802885, 
0.02290076, 0, 0, 0.03843256, 0.05144695, 0.06666667, 0.9555556, 
0.9387991, 0.9615385, 0.9310345, 0.9219219, 0.9404117, 0.953125, 
1, 0.9326923, 0.9160305, 0.9615385, 1, 0.9012811, 0.85209, 0.8888889
), .Dim = c(15L, 3L), .Dimnames = list(NULL, c("Series 1", "Series 2", 
"Series 3")), .Tsp = c(1, 15, 1), class = c("mts", "ts", "matrix"
))

You can check whether this worked by typing

obs

Now, you have a few zeros in there, which will be a problem once you take logarithms. A simple solution is to set everything that is less than a small $\epsilon$ to that $\epsilon$:

epsilon <- 0.0001
obs[obs<epsilon] <- epsilon

Now the modified rows do not sum to 1 any more. We can rectify that (although I think this might make the forecast worse):

obs <- obs/matrix(rowSums(obs),nrow=nrow(obs),ncol=ncol(obs),byrow=FALSE)

Now we transform the data as per page 35 of the presentation:

zz <- log(obs[,-ncol(obs)]/obs[,ncol(obs)])
colnames(zz) <- head(colnames(obs),-1)
zz

Load the forecast package and set a horizon of 5 time points:

library(forecast)
horizon <- 5

Now model and forecast the transformed data column by column. Here I am simply calling ets(), which will attempt to fit a state space exponential smoothing model. It turns out that it uses single exponential smoothing for all three series, but especially if you have more than 15 time periods, it may select trend models. Or if you have monthly data, explain to R that you have a potential seasonality, by using ts() with frequency=12 - then ets() will look at seasonal models.

baz <- apply(zz,2,function(xx)forecast(ets(xx),horizon=horizon)["mean"])
forecasts.transformed <- cbind(baz[[1]]$mean,baz[[2]]$mean)

Next we backtransform the forecasts as per page 38 of the presentation:

forecasts <- cbind(exp(forecasts.transformed),1)/(1+rowSums(exp(forecasts.transformed)))

Finally, let's plot histories and forecasts:

plot(obs[,1],ylim=c(0,1),xlim=c(1,nrow(obs)+horizon),type="n",ylab="")
for ( ii in 1:ncol(obs) ) {
    lines(obs[,ii],type="o",pch=19,col=ii)
    lines(forecasts[,ii],type="o",pch=21,col=ii,lty=2)
}
legend("left",inset=.01,lwd=1,col=1:ncol(obs),pch=19,legend=colnames(obs))

compositional forecasts

EDIT: a paper on compositional time series forecasting just appeared. I haven't read it, but it may be of interest.

Best Answer

Related Solutions

Solved – How to decompose a time series with multiple seasonal components

R – Solving Problems with Time Series Prediction

Related Question