当前位置:文档之家› R语言时间序列函数整理_光环大数据培训

R语言时间序列函数整理_光环大数据培训

R语言时间序列函数整理_光环大数据培训【包】library(zoo) #时间格式预处理library(xts) #同上library(timeSeires) #同上library(urca) #进行单位根检验library(tseries) #arma模型library(fUnitRoots) #进行单位根检验library(FinTS) #调用其中的自回归检验函数library(fGarch) #GARCH模型library(nlme) #调用其中的gls函数library(fArma) #进行拟合和检验【基本函数】数学函数abs,sqrt:绝对值,平方根 log, log10, log2 , exp:对数与指数函数 sin,cos,tan,asin,acos,atan,atan2:三角函数 sinh,cosh,tanh,asinh,acosh,atanh:双曲函数简单统计量sum, mean, var, sd, min, max, range, median, IQR(四分位间距)等为统计量,sort,order,rank与排序有关,其它还有ave,fivenum,mad,quantile,stem等。

#具体说明见文档1#转成时间序列类型x = rnorm(2)charvec = c(“2010-01-01”,”2010-02-01”)zoo(x,as.Date(charvec)) #包zooxts(x, as.Date(charvec)) #包xtstimeSeries(x,as.Date(charvec)) #包timeSeries#规则的时间序列,数据在规定的时间间隔内出现tm = ts(x,start = c(2010,1), frequency=12 ) #12为按月份,4为按季度,1为按年度zm = zooreg(x,start = c(2010,1), frequency=12 ) #包zooxm = as.xts(tm) #包xtssm = as.timeSeries(tm) #包timeSeries#判断是否为规则时间序列is.regular(x)#排序zoo()和xts()会强制变换为正序(按照时间名称)timeSeries不会强制排序;其结果可以根据sort函数排序,也可以采用rev()函数进行逆序;参数recordIDs,可以给每个元素(行)标记一个ID,从而可以找回原来的顺序#预设的时间有重复的时间点时xts按照升序排列timeSeries把重复部分放置在尾部;#行合并和列合并#都是按照列名进行合并,列名不同的部分用NA代替cbind()rbind()merge() 列合并#取子集xts()默认将向量做成了矩阵;其他与常规向量或者矩阵没有差别#缺失值处理na.omit(x)x[is.na(x)] = 0x[is.na(x)] = mean(x,na.rm=TRUE)x[is.na(x)] = median(x,na.rm=TRUE)na.approx(x) #对缺失值进行线性插值na.spline(x) #对缺失值进行样条插值na.locf(x) #末次观测值结转法na.trim(x, sides=”left” ) #去掉最后一个缺失值na.omit(x, “ir” ) #去掉首末位置的缺失值na.omit(x, “iz” ) #用替换首末位置的缺失值na.omit(x, “ie” ) #对首末位置的缺失值进行插值na.omit(x, method=“ie”, interp= c(“before”,”linear”,”after”) ) #可以选择插值方法,before末次观测值法,after下次观测结转法as.contiguous(x) #返回x中最长的连续无缺失值的序列片段,如果有两个等长的序列片段,则返回第一个。

#时间序列数据的显示#zoo和xts都只能按照原来的格式显示,timeSeries可以设置显示格式print(x, format= “%m/%d/%y %H:%M”) #%m表示月,%d表示天,%y表示年,%H 表示时,%M表示分钟,%A表示星期,%j表示天的序号#timeSeries也可以按照ts的格式显示print(x, style=”ts”)print(x, style=”ts”, by=”quarter”)【图形展示】plot.zoo(x)plot.xts(x)plot.zoo(x, plot.type=”single”) #支持多个时间序列数据在一个图中展示plot(x, plot.type=”single”) #支持多个时间序列数据在一个图中展示,仅对xts不行【基本统计运算】1、自相关系数、偏自相关系数等例题2.1d=scan(“sha.csv”)sha=ts(d,start=1964,freq=1)plot.ts(sha) #绘制时序图acf(sha,22) #绘制自相关图,滞后期数22pacf(sha,22) #绘制偏自相关图,滞后期数22corr=acf(sha,22) #保存相关系数cov=acf(sha,22,type = “covariance”) #保存协方差2、同时绘制两组数据的时序图d=read.csv(“double.csv”,header=F)double=ts(d,start=1964,freq=1)plot(double, plot.type = “multiple”) #两组数据两个图plot(double, plot.type = “single”) #两组数据一个图plot(double, plot.type = “single”,col=c(“red”,”green”),lty=c(1,2)) #设置每组数据图的颜色、曲线类型)3、纯随机性检验例题2.3续d=scan(“temp.csv”)Box.test(temp, type=”Ljung-Box”,lag=6)4、差分运算和滞后运算difflag5、模拟ARIMA模型的结果arima.sim(n = 100, list(ar = 0.8))plot.ts(arima.sim(n = 100, list(ar = 0.8))) #会随机产生一个包含100个随机数的时序图plot.ts(arima.sim(n = 100, list(ar = -1.1))) #非平稳,无法得到时序图。

plot.ts(arima.sim(n = 100, list(ar = c(1,-0.5))))plot.ts(arima.sim(n = 100, list(ar = c(1,0.5))))arima.sim(n = 1000, list(ar = 0.5, ma = -0.8))acf(arima.sim(n = 1000, list(ar = 0.5, ma = -0.8)),20)pacf(arima.sim(n = 1000, list(ar = 0.5, ma = -0.8)),20)【单位根检验】#方法1b=ts(read.csv(“6_1.csv”,header=T))x=b[,1]y=b[,1]#方法2:单位根检验更好的函数,加了画图的功能library(fUnitRoots)urdfTest(x)#方法3:ADF检验的一个自编函数library(urca)#…ur.df.01=function(x,lags=8){#将三种ADF检验形式汇总的函数(结果和EVIEWS不一致)res=matrix(0,5,3)colnames(res)=c(“无”,”含常数项”,”含常数项和趋势项”) rownames(res)=c(“tau统计量”,”1%临界值”,”5%临界值”, “10%临界值”,”是否稳定(1/0)”)types=c(“none”,”drift”,”trend”)for(i in 1:3){x.adf=ur.df(x,type=types[i],lags=lags,selectlags=”AIC”) x.adf.1=x.adf@teststat #统计量x.adf.2=x.adf@cval #临界值res[1,i] =x.adf.1[1]res[2:4,i]=x.adf.2[1,]res[5,i]=if( abs(res[1,i]) > abs(res[3,i]) ) 1 else 0}}#…ur.df.01(x) #对原序列进行判断【一般的ARIMA模型】d=scan(“a1.5.txt”) #导入数据prop=ts(d,start=1950,freq=1) #转化为时间序列数据plot(prop) #作时序图acf(prop,12) #作自相关图,拖尾pacf(prop,12) #作偏自相关图,1阶截尾Box.test(prop, type=”Ljung-Box”,lag=6)#纯随机性检验,p值小于5%,序列为非白噪声Box.test(prop, type=”Ljung-Box”,lag=12)( m1=arima(prop, order = c(1,0,0),method=”ML”) ) #用AR(1)模型拟合,如参数method=”CSS”,估计方法为条件最小二乘法,用条件最小二乘法时,不显示AIC。

( m2=arima(prop, order = c(1,0,0),method=”ML”, include.mean = F) ) #用AR(1)模型拟合,不含截距项。

tsdiag(m1) #对估计进行诊断,判断残差是否为白噪声summary(m1)r=m1$residuals #用r来保存残差Box.test(r,type=”Ljung-Box”,lag=6, fitdf=1)#对残差进行纯随机性检验,fitdf表示残差减少的自由度prop.fore = predict(m1, n.ahead =5) #将未来5期预测值保存在prop.fore 变量中U = prop.fore$pred + 1.96* prop.fore$se #会自动产生方差L = prop.fore$pred – 1.96* prop.fore$se #算出95%置信区间ts.plot(prop, prop.fore$pred, col=1:2) #作时序图,含预测。

lines(U, col=”blue”, lty=”dashed”)lines(L, col=”blue”, lty=”dashed”)#在时序图中作出95%置信区间——说明:运行命令arima(prop, order = c(1,0,0),method=”ML”)之后,显示:Call:arima(x = prop, order = c(1, 0, 0), method = “ML”)Coefficients:ar1 intercept0.6914 81.5509s.e. 0.0989 1.7453sigma^2 estimated as 15.51: log likelihood = -137.02, aic = 280.05注意:intercept下面的81.5509是均值,而不是截距!虽然intercept是截距的意思,这里如果用mean会更好。

相关主题