Paul Gilbert
2016-Jan-24 15:45 UTC
[R] Estimating MA parameters through arima or through package "dlm"
(Sorry for the delay in responding to this.) On 01/05/2016 06:00 AM, r-help-request at r-project.org wrote:> Date: Mon, 4 Jan 2016 11:31:22 -0500 > From: Mark Leeds<markleeds2 at gmail.com> > To: Stefano Sofia<stefano.sofia at regione.marche.it> > Cc:"r-help at r-project.org" <r-help at r-project.org> > Subject: Re: [R] Estimating MA parameters through arima or through > package "dlm" > Message-ID: > <CAHz+bWYXanjp4ksu2SzMbFRLaKH8_ASwN07s_3eKYjeyWuAG1A at mail.gmail.com> > Content-Type: text/plain; charset="UTF-8" > > Hi: I don't have time to look at the details of what you're doing but the > "equivalence" > between state space and arima ( as paul gilbert pointed out a few weeks ago > ) is not a true equivalence.To state this a bit more precisely, there is a true equivalence between the input-output equivalence classes of (linear, time invariant) state-space models with state dimension n and the input-output equivalence classes of (linear, time invariant) ARMA models with McMillan degree n. (In fact, the quotient spaces are diffeomorphic not just isomorphic.) This means you should be able to get exactly comparable results for anything that is an equivalence class invariant, including residuals and anything calculated from residuals, such as likelihood. Model roots and thus stability are also invariants, but you probably will not get comparable results for most other things involving parameters. This is not just a "statistical equivalence" as is sometimes suggested, it is an algebraic equivalence between quotient spaces. However, if you estimate a state-space model and estimate an ARMA model, there are several other things that come into play related to estimation. Comparing the two estimated models you are unlikely to find comparable results. (Typically ARMA estimation is more robust at finding the best in my experience.) Even with simulation testing of estimation starting with "true" models it is problematic to get estimated models that are equivalent. If you really want to see the equivalence you need to do a conversion of a model from one form to the other. I cannot speak to dlm, but dse was built for studying this equivalence and the users' guide has examples. If you are really interested in this topic, I recommend section 5 of http://www.bankofcanada.ca/1993/03/working-paper-199/, I realize this is getting a bit old, but if you find a more up-to-date summary I would like to hear about it. That paper also has a demonstration that the equivalent models give results that are comparable within numerical precision of computers, not just to some statistical significance.> > if you are in an area of the parameter space that the state space > formulation > can't reach, then you won't get the same parameter estimates. so, what > you're doing > might be okay or might not be, depending on whether the state space > formulation > can reach that area of the parameter space."Can't reach" is an estimation problem. This is typically a more serious problem with state-space models. The quotient spaces are diffeomorphic so in theory it should be possible to reach the same solution if you properly account for the fact that you are estimating on a smooth manifold and not a vector space. In practice you have also to worry about twisting of the parameter space and finite time for estimation, and gradients that may converge toward zero near boundaries of the manifold's charts.> there's another state space > formulation that is truly equivalent which is called the SSOE formulation > or innovations representation but > I don't know if you want to get into that. google "SSOE state space" if > you're interested.The quotient space of input-output equivalence classes for innovations form models is equivalent to the quotient space of input-output equivalence classes for non-innovations form models. You need more information to identify a non-innovations form, typically some physical understanding of the system. On the bases of only input-output data, and with no additional understanding of the physical system, there would be no reason to choose a non-innovations form for estimation. There is more discussion of this in the above mentioned summary and in the dse user's guide. BTW, estimation problems tend to be much more severe with multivariate series than with univariate series. This is not just because of the usual issues. Especially in state-space representations, the twisting of the parameter space seems to be especially bad. Paul> > > Mark > > > On Mon, Jan 4, 2016 at 9:25 AM, Stefano Sofia < > stefano.sofia at regione.marche.it> wrote: > >> >Dear list users, >> >I want to use apply a MA(2) process (x=beta1*epsilon_(t-1) + >> >beta2*epsilon_(t-1) + epsilon_(t)) to a given time series (x), and I want >> >to estimate the two parameters beta1, beta2 and the variance of the random >> >variable epsilon_(t). >> > >> >If I use >> >MA2_1 <- Arima(x, order=c(0,0,2)) >> >I get the following result >> > >> >[1] "MA2_1" >> >Series: x >> >ARIMA(0,0,2) with non-zero mean >> > >> >Coefficients: >> > ma1 ma2 intercept >> > -0.0279 0.0783 5.3737 >> >s.e. 0.0667 0.0622 0.0245 >> > >> >sigma^2 estimated as 0.1284: log likelihood=-92.63 >> >AIC=193.25 AICc=193.43 BIC=207.11 >> >[1] 0 2 0 0 1 0 0 >> > >> > From this straightforward analysis V[epsilon]=0.1284, beta1=-0.0279 and >> >beta2=0.0783. >> > >> >I also tried to use a DLM representation of ARIMA models and estimate the >> >unknown parameters by maximum likelihood through the dlm package (in >> >particular applying the example at section 3.2.6, page 115, of "Dynamic >> >Linear Models with R" by Petris, Petrone and Campagnoli: >> > >> >arma_parameters <- function(x) >> >{ >> > buildGap <- function(u) >> > { >> > gap <- dlmModARMA(ma = u[2 : 3], sigma2 = u[1]) >> > return(gap) >> > } >> > init <- c(0.005, 0.004, 0.003) >> > outMLE <- dlmMLE(x, init, buildGap) >> > dlmGap <- buildGap(outMLE$par) >> >} >> > >> >and this gives: >> >[1] "outMLE" >> >$par >> >[1] 1.00816794 0.02349296 0.02364788 >> > >> >$value >> >[1] 3089.196 >> > >> >$counts >> >function gradient >> > 10 10 >> > >> >$convergence >> >[1] 0 >> > >> >$message >> >[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" >> > >> >[1] "dlmGap" >> >$FF >> > [,1] [,2] [,3] >> >[1,] 1 0 0 >> > >> >$V >> > [,1] >> >[1,] 0 >> > >> >$GG >> > [,1] [,2] [,3] >> >[1,] 0 1 0 >> >[2,] 0 0 1 >> >[3,] 0 0 0 >> > >> >$W >> > [,1] [,2] [,3] >> >[1,] 1.00816794 0.0236848488 0.0238410337 >> >[2,] 0.02368485 0.0005564272 0.0005600964 >> >[3,] 0.02384103 0.0005600964 0.0005637899 >> > >> >$m0 >> >[1] 0 0 0 >> > >> >$C0 >> > [,1] [,2] [,3] >> >[1,] 1e+07 0e+00 0e+00 >> >[2,] 0e+00 1e+07 0e+00 >> >[3,] 0e+00 0e+00 1e+07 >> > >> >In this case >> >V[epsilon]=W[1,1]=1.00816794 >> >beta1=W[2,1]/W[1,1]=0.02349296 >> >beta2=W[3,1]/W[1,1]=0.02364788 >> > >> >I presume that these two approaches should give comparable results, but >> >this does not happen. >> >Is the model that I used correct? And does it make sense to perform this >> >kind of comparison? >> > >> >This is the log of a rainfall time series (which has already been >> >deseasonalised): >> >[1] 6.014937 4.978801 5.654592 5.616771 5.612398 5.837147 5.121580 5.832176 >> >[9] 5.205654 5.355642 5.405376 6.257859 5.516247 5.500850 4.708629 5.482304 >> >[17] 5.689684 5.727824 4.779123 5.289277 5.217107 5.976351 4.630838 >> >5.683240 >> >[25] 5.345678 5.906179 5.605434 5.497578 5.898801 5.660875 5.111988 >> >5.571013 >> >[33] 5.949340 5.374352 4.841033 5.995706 5.661223 5.458734 4.454347 >> >5.795754 >> >[41] 5.995706 5.596939 5.399971 5.908898 5.282696 5.438514 5.528635 >> >6.022721 >> >[49] 5.524257 5.519459 4.957235 5.547518 5.080783 5.411200 5.056883 >> >5.798183 >> >[57] 5.086361 5.536547 5.220356 5.141664 5.847017 5.052417 5.734635 >> >5.340419 >> >[65] 5.724238 5.634432 5.685958 5.307773 5.817706 5.134032 4.987708 >> >5.110179 >> >[73] 5.423628 5.347108 4.859037 5.556828 5.487283 5.661223 5.732370 >> >5.469325 >> >[81] 5.726848 5.419207 5.172187 5.608006 5.130490 5.586874 5.171052 >> >5.683240 >> >[89] 4.674696 5.286245 5.342813 5.370638 5.432411 5.748118 6.355239 >> >5.557986 >> >[97] 5.399067 5.222516 5.279644 5.425390 5.540871 5.917818 5.132853 >> >5.689007 >> >[105] 5.900993 5.007296 5.102911 5.778271 5.318120 5.927726 5.066385 >> >5.716699 >> >[113] 5.511815 4.714921 5.383577 5.319100 5.269403 5.354698 5.145749 >> >5.204556 >> >[121] 5.878296 5.070161 5.441552 5.213304 5.450180 5.695750 4.893352 >> >5.425390 >> >[129] 5.682559 5.487283 4.213608 5.751620 5.432411 5.379436 5.700444 >> >5.580484 >> >[137] 5.357529 5.319100 4.532599 5.603225 5.208393 5.254888 5.017280 >> >5.349961 >> >[145] 4.374498 5.187944 5.585374 5.716370 3.561046 5.119789 5.163070 >> >5.422745 >> >[153] 5.863915 5.651436 4.762174 5.655642 4.797442 5.735927 4.911183 >> >5.240688 >> >[161] 5.148076 5.477300 4.572647 5.493473 5.437644 4.854371 4.908233 >> >4.755313 >> >[169] 5.582744 5.527841 5.613128 5.211124 5.275049 5.462984 5.016617 >> >5.981919 >> >[177] 5.566817 5.094364 5.314191 5.712742 5.299317 5.452325 4.691348 >> >5.851628 >> >[185] 5.410753 5.488938 5.660179 5.900993 5.380819 5.256453 4.781641 >> >5.531807 >> >[193] 5.497578 5.274537 4.325456 5.271973 5.077047 5.258536 5.280662 >> >5.247024 >> >[201] 5.995208 4.700480 4.991113 5.457029 5.194622 5.487283 5.197391 >> >5.747161 >> >[209] 5.842094 5.372497 5.306781 5.641907 5.565286 5.259057 5.241218 >> >4.759607 >> >[217] 4.550714 5.230574 4.470495 5.664348 4.846547 5.771130 4.823502 >> >5.598422 >> >[225] 5.627621 5.547518 5.596939 5.468482 5.536940 5.606170 5.281680 >> >5.656691 >> >[233] 5.283204 5.752255 5.192401 4.550714 >> > >> > >> >Thank you for your attention and your help >> >Stefano >> > >> >