时间序列 R语言考试基本代码.doc

_典型例题HW25HW33，4HW4EXAM14,15VARlibrary(quantmod) getSymbols(AAPL,from=2008-01-03,to=2015-01-28)dim(AAPL) head(AAPL)tail(AAPL)chartSeries(AAPL,theme=white) % Obtain time plot of closing price and trading volumeda - read.csv(d-vix0411.csv,header=T)%有时候不用header参数，看head(da)的数据形式da - read.table(d-vix0411.txt,header=T)AAPL.rtn=diff(log(AAPL$AAPL.Adjusted) % Compute log returnsunrate - as.numeric(UNRATE,1) % Use a regular vector, instead of an“xts”objectlibrary(fBasics)ts.plot(ibm,main=Monthly IBM simple returns: 1968-2015) % Time plotbasicStats(ibm)apply(rtn,2,basicStats) # This command says apply basicStats to every columns in rtnd4=density(ibm)plot(d4$x,d4$y,type=l,xlab=rtn,ylab=AAPL) #density图t.test(lnIBM) % Test mean=0 vs mean .not. Zero也可以用basicStats命令去计算normalTest(lnIBM,method=jb) %testing normality of financial return seriess3=skewness(lnIBM); T - length(lnIBM)tst - s3/sqrt(6/T) % test skewnesspv - 2*pnorm(tst) %calculate p valuek4 - kurtosis(lnIBM)tst - k4/sqrt(24/T) % test excess kurtosismu - mean(sbux); v1 - var(sbux) %prediction 或者可以用t-test的信息lcl - mu-1.96*sqrt(v1)ucl - mu+1.96*sqrt(v1) c(lcl,ucl)correlationcor(sp,ibm)1 0.5785249cor(sp,ibm,method=kendall) %random copy 1 0.4172056cor(sp,ibm,method=spearman) %rank correlation1 0.58267cor(rank(ibm),rank(sp)1 0.58267x1, %show the first row of the dataX,4:5 %show the 4and5 columns of the datay=ts(x,3,frequency=252,start=c(2004,1) = Create a time-series object in Rplot(y,type=l,xlab=year,ylab=rtn)par(mfcol=c(2,1) = To put two plots on a single pagehist(y,nclass=50)%直方图tdx=(c(1:615)+11)/12+1959 %创建time plot 的横坐标plot(tdx,xt,xlab=year,ylab=temp,type=l)plot(tdx-1,zt,xlab=year,ylab=diff(temp),type=l) #注意difference使得横坐标年份变化int=cbind(x300:914,4,y,4) = Line up the two TB ratesacf(ibm)m1 - acf(ibm)names(m1)m1$acfm2 - pacf(ibm)Box.test(ibm,lag=10) % Box-Pierce Q(m) test，序列是否自相关Box.test(ibm,lag=10,type=Ljung) % Ljung-Box Q(m) testm1=lm(r3r1) % Fit a regression model with likelihood methodl AR,MA,ARMA,ARIMA modelm1=ar(x,method=mle) % Automatic AR fitting using AIC criterionnames(m1)m1$orderplot(m1$resid,type=l) 或者acf(m1$residuals) % Plot residuals of the fitted model Box.test(m1$resid,lag=10,type=Ljung) % Model checkingm2=arima(x,order=c(3,0,0) % Another approach with order giventsdiag(m2) % obtain 3 plots of model checkingtsdiag(m2,gof=24) %Model checking with 24 lags of residuals ACFc1 - c(NA,NA,NA,0,0,NA,0,0,0,NA,NA,NA) %Fixing coefficients to zero.根据t值将不显著的参数设置为0，注意一旦看出t值不显著，就要做这一步m3 - arima(unrate,order=c(11,0,0),fixed=c1)jan - rep(c(1,rep(0,11),54) %定义虚拟变量联系下面handle outlierstsdiag(m3,gof=24)source（backtest.R）backtest(m2,lhp,200,1)%比较两个model哪个更好backtest(n5,dB,orig=550,xre=X,include.mean=F)Dealing with the constant term. If there is any differencing, no constant is used.The subcommand include.mean=F in the arima command.有include.mean=F的时候估计的模型里面没有常数项% Further analysis of the fitted model.p1=c(1,-m2$coef1:3) roots=polyroot(p1)roots(大于一)Mod(roots)predict(m2,8) % Prediction 1-step to 8-step ahead and intervalsucl - p3$pred+1.96*p3$se %intervallcl - p3$pred-1.96*p3$secint which.min(m3$residuals) # Locating the largest outlier1 23 I23 I2323=1 c2 m4 tsdiag(m4,gof=24)library(fUnitRoots) %unit-root testadfTest(gdp,lag=4,type=c(c) #用ar(diff(gdp),method=”mle”)命令去获取laglibrary(forecast)Auto.arima(x) #用此命令之后要看acf图改进模型，再根据t值改进如果说auto出来的有with zero mean，一定要加上include.mean=F#Seasonal Time Series x=ts(scan(q-earn-jnj.txt),frequency=4,start=c(1960,1) % create a time series object plot(x) % Plot data with calendar time y=log(x) % Natural log transformation plot(y) % plot data 看是否有季节因素影响趋势 c1=paste(c(1:4) points(y,pch=c1) % put circles on data points. acf(y,lag.max=16) y1=as.vector(y) % Creates a sequence of data in R acf(y1,lag.max=16) dy1=diff(y1) % regular difference acf(dy1,lag.max=16) sdy1=diff(dy1,4) % seasonal difference acf(sdy1,lag.max=12) m1=arima(y1,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4) % Airline model in R. m1 tsdiag(m1) % Model checking f1=predict(m1,8) % predictionl Volatility modelsARCH,GARCH,IGARCH,GARCH-M,EGARCH,TGARCH,CHARMA,RCA,SV,LMSVlibrary(fGarch)#several types of innovational distributions cond.dist=?默认gaussian(norm) 常用std，sstdgeneralized error distribution (ged), skew normal distribution(snorm), skew Student-t (sstd), skew generalized error distribution(sged), and standardized inverse normal distribution (snig)Box.test(intc2,lag=10,type=Ljung) %检验ARCH effectm1=garchFit(garch(3,0),data=intc,trace=F) = trace=F reduces the amount of outputm2=garchFit(garch(3,0),data=intc,cond.dist=std,trace=F)m2=garchFit(arma(3,0)+garch(1,1),data=sp5,trace=F)summary(m1)plot(m1)13 rtn18411 -0.09549699 vol518411 0.744785Model checking的时候检验standardized residuals的ARCH effect，serial correlation；检验plot 的13EIGARCHsource(Egarch.R)IGARCHsource(Igarch.R)Igarch(sp5,include.mean=T) # mu is not zeroIgarch(intc,volcnt=T) # Also estimate omegaIGARCH的volatility计算很简单，用公式library(rugarch)spec1=ugarchspec(variance.model=list(model=iGARCH,garchOrder=c(1,1),mean.model=list(armaOrder=c(0,0) mm=ugarchfit(data=sp5,spec=spec1)spec1 m1 spec5 m5 spec6 m6 n3 - ugarchfit(spec=spec4,data=vw)ugarchforecast(m3,n.ahead=5) #predictionGARCH-Msource(garchM.R)m1=garchM(sp5)APGARCHTGARCH or GJR m1=garchFit(aparch(1,1),data=ibm,trace=F,delta=2,include.delta=F)m1=garchFit(aparch(1,1),data=ibm,trace=F,delta=2,include.delta=F,cond.dist=std)m1- garchFit(garch(1,1),data=ibm, trace=F, leverage=T)plot(m1)l Multivariate VARlibrary(urca)help(ca.jo) # Johansens co-integration testm2=ca.jo(x,K=2) #Test type: maximal eigenvalue statistic (lambda m