TIMESERIES ANALYSIS USING R

1、time series analysis using r by steve raper and chris chatfieldthis document is a free appendix to the 6th edition ofthe analysis of time seriesby chris chatfield, published in 2004 by chapman & hall/crc in the texts in statistical science series.the examples in this time-series text were mostly car

2、ried out using minitab and s-plus (see chapter 14 and appendix d). since the time the book was written, a software package called r has become the software choice of many students and staff and is now used extensively by the academic statistical community for statistical modelling, including time-se

3、ries modelling. r is a free updated version of s-plus. of course, many statistical software packages, including spss, sas and minitab, also contain time series analysis modules which allow the analyst to model time series as they see fit. microsoft excel even contains a command within its data analy

4、sis add-in to carry out exponential smoothing! however, the r language and environment arguably provide greater depth and flexibility in many situations. since r is used within a command-line interface, this may impose a steeper learning-curve for the new user, but the range of time series analysis

5、packages available in r, together with its publication-quality graphical capabilities, mean that it is increasingly the favoured package amongst serious time series analysts. we assume the reader is familiar with the basic use of r and indicate how to extend this to time-series analysis. the key ext

6、ension is that the time-series commands act on data called a ts object which are not just a set of numbers but have an order and a position in a cycle. for example, if the series is monthly, we may know that a particular observation is say the value in the 4th month of the second year.what is r?r is

7、 an open-source (i.e. free!) statistical software package, maintained by the user community themselves. it is distributed by cran (comprehensive r archive network) and is available for download for linux, macos x and windows from the cran web-site at . r uses the s language a

8、nd environment, developed at bell laboratories (now lucent technologies) by john chambers and colleagues, and much of the code written for s can run under r.resources in r for time series analysisa lot of resources for time series analysis are available to the r community including: several useful i

9、ndividual functions (such as plotting the sample autocorrelation and sample partial autocorrelation functions, fitting an arima model etc. for regularly spaced time series) included with the base r infrastructure additional packages for more extensive time series analysis, and for state-space models

10、 and spectral analysis time series datasets available directly in base r and in other time series packages books, on-line tutorials, and other on-line resources time-series functions available in the base r packageseveral commonly-used analysis tools for time series are available within the base r p

11、ackage, including (amongst others):acfproduces the sample autocorrelation function. user can specify maximum lags, or a vector of required lags. can also produce sample autocovariance function and sample partial autocorrelation function. pacfas acf above, but produces just the sample partial autocor

12、relation functionarimafits a sarima model of order (p,d,q)x(p,d,q), with period s. method can be chosen from:mlmaximum likelihoodcssminimising conditional sum of squarescss-mlusing conditional sum of squares to find starting values, then maximum likelihood to fit the modelpredictpredicts n steps ahe

13、ad, from any fitted model, including a time series fitted using the command arima (see above)arima.simsimulates an arima model of stated length of the order (p,d,q), with innovations having a stated variancetsdiagproduces 3 standard diagnostic charts for a fitted arima model: plot of residuals from

14、the model sample autocorrelation function of the residuals from the fitted model ljung-box portmanteau statistic for stated maximum number of lags.spectrumproduces a spectral density using one of two methods: periodogram using fast fourier transforms, optionally smoothed with daniell smoothers to be

15、 specified autoregressive fits an ar model, and computes the spectral density of the fitted model.alternatively, the command spec.prgrm can be used.time series packages in rin addition to the functions in base r, several time series analysis packages are available for specialised models, including:z

16、ooinfrastructure for both regularly- and irregularly-spaced time seriestseriescontains many specialised time series functions e.g. garch (generalised autoregressive conditional heteroscedastic) model fittingctscontinuous-time autoregressive modelsdsedynamic systems estimation tools for multivariate,

17、 time-invariant models including state-space representationsdlmbayesian and likelihood analysis of dynamic linear modelssspirtools for the specification of formulae to define and fit state-space models datasetsthere are several datasets available with base r for time series analysis, including some

18、already featured in this text. within base r, such datasets have often already been formulated into a ts (time series) class object, and the data can be recalled simply by inputting the dataset name. for example, the box-jenkins monthly airline passenger numbers (1949-1960) can be invoked as follows

19、: airpassengers jan feb mar apr may jun jul aug sep oct nov dec1949 112 118 132 129 121 135 148 148 136 119 104 1181950 115 126 141 135 125 149 170 170 158 133 114 1401951 145 150 178 163 172 178 199 199 184 162 146 1661952 171 180 193 181 183 218 230 242 209 191 172 1941953 196 196 236 235 229 243

20、264 272 237 211 180 2011954 204 188 235 227 234 264 302 293 259 229 203 2291955 242 233 267 269 270 315 364 347 312 274 237 2781956 284 277 317 313 318 374 413 405 355 306 271 3061957 315 301 356 348 355 422 465 467 404 347 305 3361958 340 318 362 348 363 435 491 505 404 359 310 3371959 360 342 406

21、396 420 472 548 559 463 407 362 4051960 417 391 419 461 472 535 622 606 508 461 390 432as can be seen, the data has been referenced by month and by year.similarly, monthly sunspot data from 1749 to 1997 can be invoked using: sunspot.month jan feb mar apr may jun jul aug sep oct nov dec1749 58.0 62.6

22、 70.0 55.7 85.0 83.5 94.8 66.3 75.9 75.5 158.6 85.21750 73.3 75.9 89.2 88.3 90.0 100.0 85.4 103.0 91.2 65.7 63.3 75.41751 70.0 43.5 45.3 56.4 60.7 50.7 66.3 59.8 23.5 23.2 28.5 44.01752 35.0 50.0 71.0 59.3 59.7 39.6 78.4 29.3 27.1 46.6 37.6 40.0 1993 59.3 91.0 69.8 62.2 61.3 49.8 57.9 42.2 22.4 56.4

23、 35.6 48.91994 57.8 35.5 31.7 16.1 17.8 28.0 35.1 22.5 25.7 44.0 18.0 26.21995 24.2 29.9 31.1 14.0 14.5 15.6 14.5 14.3 11.8 21.1 9.0 10.01996 11.5 4.4 9.2 4.8 5.5 11.8 8.2 14.4 1.6 0.9 17.9 13.31997 5.7 7.6 8.7 15.5 18.5 12.7 10.4 24.4 51.3 22.8 39.0 41.2other regular time series can be read into r

24、and turned into a ts object using the ts command for example:ts(x, frequency = 4, start = c(1959, 2) transforms the vector x into a quarterly time series object (frequency = 4), starting in quarter 2 of the year 1959.time series packages have their own inherent datasets and time series -objects.on-l

25、ine and off-line resourcesseveral springer texts in statistics cover time series analysis using r, including:time series analysis & its applications: with r examples (3rd ed) robert shumway and david stoffer, springer (2011)(in addition stoffers own web-site includes a useful r time series tutorial

26、at /stoffer/tsa2/r_time_series_quick_fix.htm)introductory time series with r examples paul cowpertwait and andrew metcalfe, springer (2009)time series analysis in r examples (2nd ed) jonathan cryer and kung-sik chan, springer (2009)three examplesthree examples have been select

27、ed to illustrate the use of r in time series analysis and to provide further guidance for the reader. they correspond to examples 14.1, 14.2, and 14.3 in the book. this appendix concentrates on the r commands. for further discussion of the modelling process, see the book.example 1.monthly air temper

28、ature at recifetable 14.1 in the book shows the air temperature at recife in brazil in successive months over a 10-year period. the objective of our analysis is to describe and understand the data.the first step is to import the data into r and to produce an object of class ts:recife=ts(read.csv(e:r

29、ecife.csv,header=false),start=1953,frequency=12)recife jan feb mar apr may jun jul aug sep oct nov dec1953 26.8 27.2 27.1 26.3 25.4 23.9 23.8 23.6 25.3 25.8 26.4 26.91954 27.1 27.5 27.4 26.4 24.8 24.3 23.4 23.4 24.6 25.4 25.8 26.71955 26.9 26.3 25.7 25.7 24.8 24.0 23.4 23.5 24.8 25.6 26.2 26.51956 2

30、6.8 26.9 26.7 26.1 26.2 24.7 23.9 23.7 24.7 25.8 26.1 26.51957 26.3 27.1 26.2 25.7 25.5 24.9 24.2 24.6 25.5 25.9 26.4 26.91958 27.1 27.1 27.4 26.4 25.5 24.7 24.3 24.4 24.8 26.2 26.3 27.01959 27.1 27.5 26.2 28.2 27.1 25.4 25.6 24.5 24.7 26.0 26.5 26.81960 26.3 26.7 26.6 25.8 25.2 25.1 23.3 23.8 25.2

31、25.5 26.4 26.71961 27.0 27.4 27.0 26.3 25.9 24.6 24.1 24.3 25.2 26.3 26.4 26.7table 14.1:we now plot the data, as seen below:plot(recife,ylab=temperature (degree c), xlab=year,main=recife, brazil temperature data)the plot exhibits regular seasonal variation with little or no trend, as we would expec

32、t a priori. the correlogram is produced using the r command: acf(ts(recife,freq=1),lag.max=40,main=autocorrelation function for recife data,ylim=c(-1,1)(note that the r function acfs default plot gives an x-axis in terms of the frequency (in this case 12) of the time series object recife. the autoco

33、rrelation function for lag 12 months is then labelled 1 year. re-stating the frequency as 1 allows a plot of the autocorrelation function as we would expect it)the correlogram identifies the obvious seasonal variation, with high positive autocorrelations at lags 12, 24, we can remove the seasonality

34、 in the data by calculating monthly averages (see table below) and subtracting them from the raw data:monthav. temp. 1953-161 (oc)january26.82february27.08march26.70april26.32may25.60june24.62july24.00august23.98september24.98october25.83november26.28december26.74the resulting, deseasonalised data i

35、s shown below:with the above (plot of the) sample ac.f (ignore the word periodogram in the title)the raw periodogram for the deseasonalised data is given below followed by three examples of smoothed periodograms using different daniell filters (3,3), (5,5), and (7,7) to illustrate the smoothing effe

36、ct:example 2yield on short-term government securitieswe first import and plot the data:yield1-read.csv(f:ma30085 50085 time serieschatfield data filesyield.csv,header=false)yield yield jan feb mar apr may jun jul aug sep oct nov dec1950 2.22 2.23 2.22 2.20 2.09 1.97 2.03 1.98 1.94 1.79 1.74 1.861951

37、 1.78 1.72 1.79 1.82 1.89 1.99 1.89 1.83 1.71 1.70 1.97 2.211952 2.36 2.41 2.92 3.15 3.26 3.51 3.48 3.16 3.01 2.97 2.88 2.911953 3.45 3.29 3.17 3.09 3.02 2.99 2.97 2.94 2.84 2.85 2.86 2.891954 2.93 2.93 2.87 2.82 2.63 2.33 2.22 2.15 2.28 2.28 2.06 2.541955 2.29 2.66 3.03 3.17 3.83 3.99 4.11 4.51 4.6

38、6 4.37 4.45 4.581956 4.58 4.76 4.89 4.65 4.51 4.65 4.52 4.52 4.57 4.65 4.74 5.101957 5.00 4.74 4.79 4.83 4.80 4.83 4.77 4.80 5.38 6.18 6.02 5.911958 5.66 5.42 5.06 4.70 4.73 4.64 4.62 4.48 4.43 4.33 4.32 4.301959 4.26 4.02 4.06 4.08 4.09 4.14 4.15 4.20 4.30 4.26 4.15 4.271960 4.69 4.72 4.92 5.10 5.2

39、0 5.56 6.08 6.13 6.09 5.99 5.58 5.591961 5.42 5.30 5.44 5.32 5.21 5.47 5.96 6.50 6.48 6.00 5.83 5.911962 5.98 5.91 5.64 5.49 5.43 5.33 5.22 5.03 4.74 4.55 4.68 4.531963 4.67 4.81 4.98 5.00 4.94 4.84 4.76 4.67 4.51 4.42 4.53 4.701964 4.75 4.90 5.06 4.99 4.96 5.03 5.22 5.47 5.45 5.48 5.57 6.331965 6.6

40、7 6.52 6.60 6.78 6.79 6.83 6.91 6.93 6.65 6.53 6.50 6.691966 6.58 6.42 6.79 6.82 6.76 6.88 7.22 7.41 7.27 7.03 7.09 7.181967 6.69 6.50 6.46 6.35 6.31 6.41 6.60 6.57 6.59 6.80 7.16 7.511968 7.52 7.40 7.48 7.42 7.53 7.75 7.80 7.63 7.51 7.49 7.64 7.921969 8.10 8.18 8.52 8.56 9.00 9.34 9.04 9.08 9.14 8.

41、99 8.96 8.861970 8.79 8.62 8.29 8.05 8.00 7.89 7.48 7.31 7.42 7.51 7.71 7.99the sample ac.f is produced using the command:yieldacf-acf(yield1,lag.max=50,main=correlogram for monthly yield from british short-term government securities,ylim=c(-1,1)yieldacfvalues of the ac.f are listed below, because t

42、hey differ from the values quoted in chatfield 6th edition, namelylagchatfieldthis analysis10.970.9920.940.9730.920.9540.890.9350.850.91240.340.52: 1, 1.0000000 2, 0.9855462 3, 0.9680912 4, 0.9506049 5, 0.9317940 6, 0.9126326 7, 0.8922248 8, 0.8707009 9, 0.849594610, 0.829579211, 0.809158512, 0.7883

43、94513, 0.767441814, 0.744867815, 0.722140816, 0.700751117, 0.679276218, 0.658558219, 0.639071920, 0.617811721, 0.596046422, 0.574864523, 0.554477224, 0.536623625, 0.519586426, 0.504567827, 0.492092128, 0.483704029, 0.477261530, 0.471433531, 0.466190832, 0.460994633, 0.455473534, 0.450016935, 0.44401

44、4736, 0.437533337, 0.4302790these differences are being investigated by the author and are thought to be caused by an error in the data as originally listed. the correlogram is plotted below and shows the slow decay typical of a non-stationary series.we now take first differences of the data and plo

45、t the differenced series:diffyield=diff(yield, lag=1)plot(diffyield,main=differenced monthly yield from british short-term government securities,ylab=percent,xlab=year,ylim=c(-1,1)the sample ac.f of the first differences is produced via:diffyieldacf diffyieldacf$acf, , 1 1, 1.000000e+00 2, 3.424154e

46、-01 3, 4.511240e-02 4, 1.285514e-01 5, 2.119021e-02 6, 1.943120e-02 7, -1.224694e-02 8, -9.657781e-02 9, -1.381191e-0110, -8.739276e-0211, -9.216505e-0212, -1.436796e-0213, 8.897098e-0214, -4.771946e-0215, -1.274294e-0116, -4.375536e-0217, -4.926974e-0218, -8.955549e-0219, 3.081316e-0220, 1.091890e-

47、0121, 2.312124e-0222, -7.968747e-0223, -7.428385e-0224, 4.312129e-0525, -1.045158e-0126, -1.503248e-0127, -2.002343e-0128, -1.303668e-0129, -7.093623e-0230, -4.409803e-0231, -9.224070e-03fitting an arima(0,1,1) model to the original data, using the css-ml and css methods:fit-arima(yield,order=c(0,1,

48、1), method=css-ml)call:arima(x = yield, order = c(0, 1, 1), method = css-ml)coefficients: ma1 0.4350s.e. 0.0652sigma2 estimated as 0.03576: log likelihood = 61.76, aic = -119.52fit-arima(yield,order=c(0,1,1), method=css)call:arima(x = yield, order = c(0, 1, 1), method = css)coefficients: ma1 0.4373s

49、.e. 0.0653sigma2 estimated as 0.03576: part log likelihood = 61.86note that the results are very similar in this case.example 3airline passenger datamonthly totals of international airline passengers in the us are available as a time series object (named airpassengers) within the r base package, and

50、 are tabulated in table 14.3 in the book.we start by plotting the data, using the command:plot(airpassengers,ylab=number of passengers (in thousands), xlab=year)the seasonal pattern and positive trend in the data is clear. the magnitude of the seasonal variation increases at a similar rate to the ye

51、arly mean levels, indicating a multiplicative seasonal model is appropriate.we attempt to fit an arima model to the data, using r.firstly, the multiplicative seasonality of the data suggests we transform the data by taking logarithms, in order to make the seasonality additive.logap-log(airpassengers)plot(logap,ylab=logarithm of number of passengers (in thousands), xlab=year)following the log transformation, the plot (not shown here) shows the size of the seasonal variation is now roughly constant. we can now plot the ac.f of the raw logged data, together with those for the first diff