Time series ARMA model estimation results different from textbook

I tried to estimate an ARMA(4,3) model using armax function in the system identification toolbox. I used the same data set, but the result was totally different from the textbook Pandit and Wu,1990, Time Series and System Analysis with Applications. The data I used is the sunspot data in Table A2 in the textbook (page487). I double checked to make sure the data is exactly the same as in the textbook.

The estimated ARMA(4,3) model in the textbook in page 170 is:

A(q)=1-0.42(+-1)q^-1-0.28(+-0.75)q^-2+0.08(+-0.49)q^-3+0.33(+-0.5)q^-4.
C(q)=1+0.89(+-1)q^-1+0.31(+-0.76)q^-2-0.14(+-0.23)q^-3.

However, the estimation result of armax is:

A(q)=1-1.572(+-0.3967)q^-1+0.1358(+-0.9669)q^-2+1.046(+-0.8845)q^-3-0.6024(+-0.3122)q^-4.
C(q)=1-0.3253(+-0.4095)q^-1-0.756(+-0.5163)q^-2 +0.08137(+-0.168)q^-3.

My code is as follows. My question is which result is correct?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Read and plot raw data (rawdata has 176 points in one column)
rawdata=[80.9 83.4 47.7 47.8 30.7 12.2 9.6 ...
10.2 32.4 47.6 54 62.9 85.9 61.2 45.1 36.4 ...
20.9 11.4 37.8 69.8 106.1 100.8 81.6 66.5 ...
34.8 30.6 7 19.8 92.5 154.4 125.9 84.8 ...
68.1 38.5 22.8 10.2 24.1 82.9 132 130.9 ...
118.1 89.9 66.6 60 46.9 41 21.3 16 ...
6.4 4.1 6.8 14.5 34 45 43.1 47.5 ...
42.2 28.1 10.1 8.1 2.5 0 1.4 ...
5 12.2 13.9 35.4 45.8 41.1 30.4 23.9 ...
15.7 6.6 4 1.8 8.5 16.6 36.3 49.7 ...
62.5 67 71 47.8 27.5 8.5 13.2 56.9 ...
121.5 138.3 103.2 85.8 63.2 36.8 24.2 10.7 ...
15 40.1 61.5 98.5 124.3 95.9 66.5 64.5 ...
54.2 39 20.6 6.7 4.3 22.8 54.8 93.8 ...
95.7 77.2 59.1 44 47 30.5 16.3 7.3 ...
37.3 73.9 139.1 111.2 101.7 66.3 44.7 17.1 ...
11.3 12.3 3.4 6 32.3 54.3 59.7 63.7 ...
63.5 52.2 25.4 13.1 6.8 6.3 7.1 35.6 ...
73 84.9 78 64 41.8 26.2 26.7 12.1 9.5 ...
2.7 5 24.4 42 63.5 53.8 62 48.5 43.9 ...
18.6 5.7 3.6 1.4 9.6 47.4 57.1 103.9 ...
80.6 63.6 37.6 26.1 14.2 5.8 16.7]';

meanvalue=mean(rawdata);

%Fit an ARMA model
order=[4 3];
data=iddata(rawdata-meanvalue,[],1);
model=armax(data,order);

%Show the fitted model
present(model);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Walter Roberson]

I really have no idea, but the 1.572 coefficient is highly suggestive of Pi/2 and thus hints at a possible difference in whether the values are being interpreted in some absolute unit (e.g., Hz) or in some relative unit such as "cycles per radian"

[Wayne King]

I don't have this book, but are you sure you did not leave out some important detail about how they preprocessed the data. Are you able to post a copy of the relevant pages?

Answer 1

[Wayne King]

Hi Yong, I don't know what to say, but an ARMA(4,3) model fitted in R with arima is not close to what the book lists either. Perhaps there is a problem with the book? Did they use MATLAB, R, none of the above?

Answer 2

[Yong]

Hi Walter & Wayne, thank you for the comments. I don't know how to post pictures here. A scanned copy of the pages is here: http://dl.dropbox.com/u/24398626/Question%20on%20Sunspot%20ARMA%20modeling%20using%20armax.pdf

For ARMA(2,1), ARMA(2,0), ARMA(1,0) models, armax outputs are the same as the textbook. The output for ARMA(4,3) is different.

Answer 3

[Yong]

Hi Wayne, again, thanks for the reply. I've tried two other data sets in the book. The same thing happened. The ARMA(2,1) model fitted with armax is the same as the textbook, but the ARMA(4,3) and ARMA(6,5) models fitted are not.

The textbook was published in 1983, and they used their own FORTRAN code for fitting. It was used to teach undergraduate and graduate student at the University of Wisconsin and the University of Michigan.

One more question: Is the ARMA(4,3) model fitted in R with arima the same as the model fitted in matlab with armax? If it is, I will assume that the result from the textbook is not so accurate.

[Wayne King]

Yes, the R and Matlab models are close if I set the parameters the same.

[Yong]

Thanks

Answer 4

[Gurudatha Pai]

Assuming that the there is no typo in the data (I dint check), I tried your code through system id GUI of Matlab system id toolbox. It give the same values as you have typed here (obviously, right!) for 'zero' initialization and gives (very) different answers for other initializations. So, I am not sure about the books answer! I am not familiar with this book. I am just wondering what the "initial values" are supposed to be! Is this something on any of the recursive system id algorithms which are using initial values of the parameters? I don't see that you have used these initial values anywhere. You might want to use these initial values and see what happens to the result. I don't know if system id toolbox allows you to use non-zero initial values.

If I get a chance, I will play more with your data and may be I will have a better answer later! Finally, I would take it for granted that any computation tool (Matlab, R or Fortan for that matter) should give visibly same answer (unless some numerical or theoretical exception is happening).

Answer 5

[Vibhav Gaur]

The book provides a hint as to why they may be different (at the end of the sunspot model fitting example (page 169 in my book):

"It should be noted ... that the actual results obtained by the computer programs listed at the end may be somewhat different due to a variety of reasons such as subjective choices of convergence tolerances, initial values, method of computing 's etc."

I wonder if checking their FORTRAN code at the end and making sure all settings are the same will result in estimation of similar (if not the same) parameters.