R在水文时间序列分析的应用.docx

1、R在水文时间序列分析的应用R在水文时间序列分析的应用自回归滑动平均模型Autoregressive Models - AR(p)ar stats Fit Autoregressive Models to Time SeriesDescriptionFit an autoregressive time series model to the data, by default selecting the complexity by AIC.Usagear(x, aic = TRUE, order.max = NULL, method=c(yule-walker, burg, ols, mle, y

2、w), na.action, series, .)ArgumentsxA univariate or multivariate time series.aicLogical flag. If TRUE then the Akaike Information Criterion is used to choose the order of the autoregressive model. If FALSE, the model of order order.max is fitted.order.maxMaximum order (or order) of model to fit. Defa

3、ults to the smaller of N-1 and 10*log10(N) where N is the number of observations except for method=mle where it is the minimum of this quantity and 12.methodCharacter string giving the method used to fit the model. Must be one of the strings in the default argument (the first few characters are suff

4、icient). Defaults to yule-walker.na.actionfunction to be called to handle missing values.seriesnames for the series. Defaults to deparse(substitute(x).在概率论中，一个时间序列是一串随机变量。在统计学中，这样一些变量都会受时间影响：比如每天在变的股票价格，每月一测的空气温度，每分钟病人的心率等等数据：北美五大湖之一的Lake Huron的1875-1972年每年的水位值这个时间序列大致的图像：plot(LakeHuron, ylab = , ma

5、in = Level of Lake Huron)AR（1）模型：x-LakeHuronop - par(mfrow=c(2,1)y - filter(x,.8,method=recursive)plot(y, main=AR(1), ylab=)acf( y, main = paste( p =, signif( dwtest( y 1 ) $ p.value, 3 ) )par(op)ACF和PCF图op - par(mfrow = c(3,1), mar = c(2,4,1,2)+.1)acf(x, xlab = )pacf(x, xlab = )spectrum(x, xlab = ,

6、 main = )par(op)AR(p)模型使用Yule-walker法得出估计的参数值 y 1.1319*x98-0.1319*x97 1 579.9692参考书目：Introductory Time Series with R，Analysis of Time Series Data Using R，Time Series Analysis and Its Applications - with R examples，Time Series Analysis and Its Applications - with R examples参考网站：http:/zoonek2.free.fr/

7、UNIX/48_R/15.html#2MA (Moving Average models) Here is a simple way of building a time series from a white noise: just perform a Moving Average (MA) of this noise. n - 200 x - rnorm(n)y - ( x2:n + x2:n-1 ) / 2op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=white noise)plot(ts(y), xla

8、b=, ylab=MA(1)acf(y, main=)par(op)n - 200 x - rnorm(n)y - ( x1:(n-3) + x2:(n-2) + x3:(n-1) + x4:n )/4op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=white noise)plot(ts(y), xlab=, ylab=MA(3)acf(y, main=)par(op)You can also compute the moving average with different coefficients. n -

9、200 x - rnorm(n)y - x2:n - x1:(n-1)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=white noise)plot(ts(y), xlab=, ylab=momentum(1)acf(y, main=)par(op)n - 200 x - rnorm(n)y - x3:n - 2 * x2:(n-1) + x1:(n-2)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=white noise)pl

10、ot(ts(y), xlab=, ylab=Momentum(2)acf(y, main=)par(op)Instead of computing the moving average by hand, you can use the filter function. n - 200 x - rnorm(n)y - filter(x, c(1,-2,1)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=White noise)plot(ts(y), xlab=, ylab=Momentum(2)acf(y, na.

11、action=na.pass, main=)par(op)TODO: the side=1 argument.AR (Auto-Regressive models) Another means of building a time series is to compute each term by adding noise to the preceding term: this is called a random walk. For instance, n - 200x - rep(0,n)for (i in 2:n) xi - xi-1 + rnorm(1)This can be writ

12、ten, more simply, with the cumsum function. n - 200x - rnorm(n)y - cumsum(x)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=)plot(ts(y), xlab=, ylab=AR(1)acf(y, main=)par(op)More generally, one can consider X(n+1) = a X(n) + noise.This is called an auto-regressive model, or AR(1), b

13、ecause one can estimate the coefficients by performing a regression of x against lag(x,1). n - 200a - .7x - rep(0,n)for (i in 2:n) xi - a*xi-1 + rnorm(1)y - x-1x - x-nr - lm( y x -1)plot(yx)abline(r, col=red)abline(0, .7, lty=2)More generally, an AR(q) process is a process in which each term is a li

14、near combination of the q preceding terms and a white noise (with fixed coefficients). n - 200x - rep(0,n)for (i in 4:n) xi - .3*xi-1 -.7*xi-2 + .5*xi-3 + rnorm(1)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=AR(3)acf(x, main=, xlab=)pacf(x, main=, xlab=)par(op)You can also simula

15、te those models with the arima.sim function. n - 200x - arima.sim(list(ar=c(.3,-.7,.5), n)op - par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)plot(ts(x), xlab=, ylab=AR(3)acf(x, xlab=, main=)pacf(x, xlab=, main=)par(op)PACF The partial AutoCorrelation Function (PACF) provides an estimation of the coefficients

16、of an AR(infinity) model: we have already seen it on the previous examples. It can be easily computed from the autocorrelation function with the Yule-Walker equations. Yule-Walker Equations To compute the auto-correlation function of an AR(p) process whose coefficients are known, (1 - a1 B - a2 B2 -

17、 . - ap Bp) Y = Zwe just have to compute the first autocorrelations r1, r2, ., rp, and then use the Yule-Walker equations: r(j) = a1 r(j-1) + a2 r(j-2) + . + ap r(j-p).You can also use them in the other direction to compute the coefficients of an AR process from its autocorrelations. DescriptionFit

18、an ARMA model to a univariate time series by conditional least squares. For exact maximum likelihood estimation seearima0.Usagearma(x, order = c(1, 1), lag = NULL, coef = NULL, include.intercept = TRUE, series = NULL, qr.tol = 1e-07, .)Argumentsxa numeric vector or time series.ordera two dimensional

19、 integer vector giving the orders of the model to fit.order1corresponds to the AR part andorder2to the MA part.laga list with componentsarandma. Each component is an integer vector, specifying the AR and MA lags that are included in the model. If both,orderandlag, are given, only the specification f

20、romlagis used.coefIf given this numeric vector is used as the initial estimate of the ARMA coefficients. The preliminary estimator suggested in Hannan and Rissanen (1982) is used for the default initialization.include.interceptShould the model contain an intercept?seriesname for the series. Defaults

21、 todeparse(substitute(x).qr.tolthetolargument forqrwhen computing the asymptotic standard errors ofcoef.Examplesdata(tcm) r - diff(tcm10y)summary(r.arma - arma(r, order = c(1, 0)summary(r.arma - arma(r, order = c(2, 0)summary(r.arma - arma(r, order = c(0, 1)summary(r.arma - arma(r, order = c(0, 2)su

22、mmary(r.arma - arma(r, order = c(1, 1)plot(r.arma)data(nino)s - nino3.4summary(s.arma - arma(s, order=c(20,0)summary(s.arma - arma(s, lag=list(ar=c(1,3,7,10,12,13,16,17,19),ma=NULL)acf(residuals(s.arma), na.action=na.remove)pacf(residuals(s.arma), na.action=na.remove)summary(s.arma - arma(s, lag=list(ar=c(1,3,7,10,12,13,16,17,19),ma=12)summary(s.arma - arma(s, lag=list(ar=c(1,3,7,10,12,13,16,17),ma=12)plot(s.arma) （本资料素材和资料部分来自网络，仅供参考。请预览后才下载，期待您的好评与关注！）