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

R在水文时间序列分析的应用自回归滑动平均模型AutoregressiveModels-AR(p)ar{stats}FitAutoregressiveModelstoTimeSeriesDescriptionFitanautoregressivetimeseriesmodeltothedata,bydefaultselectingthecomplexitybyAIC.Usagear(x,aic=TRUE,order.max=NULL,method=c("yule-walker","burg","ols","mle","yw"),na.action,series,...)ArgumentsxAunivariateormultivariatetimeseries.aicLogicalflag.IfTRUEthentheAkaikeInformationCriterionisusedtochoosetheorderoftheautoregressivemodel.IfFALSE,themodeloforderorder.maxisfitted.order.maxMaximumorder(ororder)ofmodeltofit.DefaultstothesmallerofN-1and10*log10(N)whereNisthenumberofobservationsexceptformethod="mle"whereitistheminimumofthisquantityand12.methodCharacterstringgivingthemethodusedtofitthemodel.Mustbeoneofthestringsinthedefaultargument(thefirstfewcharactersaresufficient).Defaultsto"yule-walker".na.actionfunctiontobecalledtohandlemissingvalues.seriesnamesfortheseries.Defaultstodeparse(substitute(x)).在概率论中，一个时间序列是一串随机变量。在统计学中，这样一些变量都会受时间影响：比如每天在变的股票价格，每月一测的空气温度，每分钟病人的心率等等数据：北美五大湖之一的LakeHuron的1875-1972年每年的水位值这个时间序列大致的图像：plot(LakeHuron,ylab="",main="LevelofLakeHuron")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="",main="")par(op)AR(p)模型使用Yule-walker法得出估计的参数值y<-ar(x,aic=TRUE,method="yule-walker")regr=ar.ols(x,order=2,demean=FALSE,intercept=FALSE)regr结果：Call:ar.ols(x=x,order.max=2,demean=FALSE,intercept=FALSE)Coefficients:121.1319-0.1319Orderselected2sigma^2estimatedas0.5281预测1973值>1.1319*x[98]-0.1319*x[97][1]579.9692参考书目：IntroductoryTimeSerieswithR，AnalysisofTimeSeriesDataUsingR，TimeSeriesAnalysisandItsApplications--withRexamples，TimeSeriesAnalysisandItsApplications--withRexamples参考网站：http://zoonek2.free.fr/UNIX/48_R/15.html#2MA(MovingAveragemodels)Hereisasimplewayofbuildingatimeseriesfromawhitenoise:justperformaMovingAverage(MA)ofthisnoise.n<-200x<-rnorm(n)y<-(x[2:n]+x[2:n-1])/2op<-par(mfrow=c(3,1),mar=c(2,4,2,2)+.1)plot(ts(x),xlab="",ylab="whitenoise")plot(ts(y),xlab="",ylab="MA(1)")acf(y,main="")par(op)n<-200x<-rnorm(n)y<-(x[1:(n-3)]+x[2:(n-2)]+x[3:(n-1)]+x[4:n])/4op<-par(mfrow=c(3,1),mar=c(2,4,2,2)+.1)plot(ts(x),xlab="",ylab="whitenoise")plot(ts(y),xlab="",ylab="MA(3)")acf(y,main="")par(op)Youcanalsocomputethemovingaveragewithdifferentcoefficients.n<-200x<-rnorm(n)y<-x[2:n]-x[1:(n-1)]op<-par(mfrow=c(3,1),mar=c(2,4,2,2)+.1)plot(ts(x),xlab="",ylab="whitenoise")plot(ts(y),xlab="",ylab="momentum(1)")acf(y,main="")par(op)n<-200x<-rnorm(n)y<-x[3:n]-2*x[2:(n-1)]+x[1:(n-2)]op<-par(mfrow=c(3,1),mar=c(2,4,2,2)+.1)plot(ts(x),xlab="",ylab="whitenoise")plot(ts(y),xlab="",ylab="Momentum(2)")acf(y,main="")par(op)Insteadofcomputingthemovingaveragebyhand,youcanusethe"filter"function.n<-200x<-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="Whitenoise")plot(ts(y),xlab="",ylab="Momentum(2)")acf(y,na.action=na.pass,main="")par(op)TODO:the"side=1"argument.AR(Auto-Regressivemodels)Anothermeansofbuildingatimeseriesistocomputeeachtermbyaddingnoisetotheprecedingterm:thisiscalledarandomwalk.Forinstance,n<-200x<-rep(0,n)for(iin2:n){x[i]<-x[i-1]+rnorm(1)}Thiscanbewritten,moresimply,withthe"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)Moregenerally,onecanconsiderX(n+1)=aX(n)+noise.Thisiscalledanauto-regressivemodel,orAR(1),becauseonecanestimatethecoefficientsbyperformingaregressionofxagainstlag(x,1).n<-200a<-.7x<-rep(0,n)for(iin2:n){x[i]<-a*x[i-1]+rnorm(1)}y<-x[-1]x<-x[-n]r<-lm(y~x-1)plot(y~x)abline(r,col='red')abline(0,.7,lty=2)Moregenerally,anAR(q)processisaprocessinwhicheachtermisalinearcombinationoftheqprecedingtermsandawhitenoise(withfixedcoefficients).n<-200x<-rep(0,n)for(iin4:n){x[i]<-.3*x[i-1]-.7*x[i-2]+.5*x[i-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)Youcanalsosimulatethosemodelswiththe"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)PACFThepartialAutoCorrelationFunction(PACF)providesanestimationofthecoefficientsofanAR(infinity)model:wehavealreadyseenitonthepreviousexamples.Itcanbeeasilycomputedfromtheautocorrelationfunctionwiththe"Yule-Walker"equations.Yule-WalkerEquationsTocomputetheauto-correlationfunctionofanAR(p)processwhosecoefficientsareknown,(1-a1B-a2B^2-...-apB^p)Y=Zwejusthavetocomputethefirstautocorrelationsr1,r2,...,rp,andthenusetheYule-Walkerequations:r(j)=a1r(j-1)+a2r(j-2)+...+apr(j-p).YoucanalsousethemintheotherdirectiontocomputethecoefficientsofanARprocessfromitsautocorrelations.DescriptionFitanARMAmodeltoaunivariatetimeseriesbyconditionalleastsquares.Forexactmaximumlikelihoodestimationsee

arima0.Usagearma(x,order=c(1,1),lag=NULL,coef=NULL,ercept=TRUE,series=NULL,qr.tol=1e-07,...)Argumentsxanumericvectorortimeseries.orderatwodimensionalintegervectorgivingtheordersofthemodeltofit.

order[1]

correspondstotheARpartand

order[2]

totheMApart.lagalistwithcomponents

ar

and

ma.Eachcomponentisanintegervector,specifyingtheARandMAlagsthatareincludedinthemodel.Ifboth,

order

andlag,aregiven,onlythespecificationfrom

lag

isused.coefIfgiventhisnumericvectorisusedastheinitialestimateoftheARMAcoefficients.ThepreliminaryestimatorsuggestedinHannanandRissanen(1982)isusedforthedefaulterceptShouldthemodelcontainanintercept?seriesnamefortheseries.Defaultsto

deparse(substitute(x)).qr.tolthe

tol

argumentfor

HYPERLINK"/Rdoc/li