R語言常用上機命令分功能整理——時間序列分析為主

1、R語言常用上機命令分功能整理一一時間序列分析為主第一講應用實例?R的基本界面是一個交互式命令窗口，命令提示符是一個大于號，命令的結果馬上顯示在命令下面。?S命令主要有兩種形式：表達式或賦值運算(用W或者表示)。在命令提示符后鍵入一個表達式表示計算此表達式并顯示結果。賦值運算把賦值號右邊的值計算出來賦給左邊的變量。?可以用向上光標鍵來找回以前運行的命令再次運行或修改后再運行。?S是區(qū)分大小寫的，所以x和X是不同的名字。我們用一些例子來看R軟件的特點。假設我們已經(jīng)進入了R的交互式窗口。如果沒有打開的圖形窗口，在R中，用：x11()可以打開一個作圖窗口。然后，輸入以下語句：x1=0:100x2=x1

2、*2*pi/100y=sin(x2)plot(x2,y,type="l")這些語句可以繪制正弦曲線圖。其中，“建賦值運算符。0:100表示一個從0到100的等差數(shù)列向量。第二個語句可以看出，我們可以對向量直接進行四則運算，計算得到的x2是向量x1的所有元素乘以常數(shù)2*pi/100的結果。從第三個語句可看到函數(shù)可以以向量為輸入，并可以輸出一個向量，結果向量y的每一個分量是自變量x2的每一個分量的正弦函數(shù)值。plot（x2,y,type="l",main="畫圖練習",sub="好好練",xlab="x軸&q

3、uot;,ylab='y軸'）有關作圖命令plot的詳細介紹可以在R中輸入help(plot)數(shù)學函數(shù)abs,sqrt:絕對值，平方根log,log10,log2,exp：對數(shù)與指數(shù)函數(shù)sin,cos,tan,asin,acos,atan,atan2：三角函數(shù)sinh,cosh,tanh,asinh,acosh,atanh：雙曲函數(shù)簡單統(tǒng)計量sum,mean,var,sd,min,max,range,median,IQR(四分位間品巨)等為統(tǒng)計量，sort,order,rank與排序有關，其它還有ave,fivenum,mad,quantile,stem等。下面我們看一看S的統(tǒng)

4、計功能：> marks<-c(10,6,4,7,8)> mean(marks)> sd(marks)> min(marks)> max(marks)第一個語句輸入若干數(shù)據(jù)到一個向量，函c()用來把數(shù)據(jù)組合為一個向量。后面用了幾個函數(shù)來計算數(shù)據(jù)的均值、標準差、最小值、最大值?？梢园讶舾尚忻畋４嬖谝粋€文本文件中，然后用source函數(shù)來運行整個文件：> source("C:/l.R")注意字符串中的反斜杠。例:計算6,4,7,8,10的均值和標準差，把若干行命令保存在一個文本文件(比如C:1.R)中，然后用source函數(shù)來運行整個文

5、件。a<-c(10,6,4,7,8)b<-mean(a)c<-sd(a)source("C:/1.R")時間序列數(shù)據(jù)的輸入使用函數(shù)tsts(1:10,frequency=4,start=c(1959,2)print(ts(1:10,frequency=7,start=c(12,2),calendar=TRUE)a<-ts(1:10,frequency=4,start=c(1959,2)plot(a)將外部數(shù)據(jù)讀入Rread.csv默認header=TRUE,也就是第一行是標簽，不是數(shù)據(jù)。read.table默認header=FALSE將R中的數(shù)據(jù)輸出

6、writewrite.tablewrite.csv第二講1 .繪制時序圖、自相關圖例題2.1d=scan("sha.csv")sha=ts(d,start=1964,freq=1)plot.ts(sha)#繪制時序圖acf(sha,22)#繪制自相關圖，滯后期數(shù)22pacf(sha,22)#繪制偏自相關圖，滯后期數(shù)22corr=acf(sha,22)#保存相關系數(shù)cov=acf(sha,22,type="covariance")琳存協(xié)方差圖的保存，單擊選中圖，在菜單欄選中“文件”，再選“另存為”。同時顯示多個圖：用x11()命令生成一個空白圖，再輸入作圖

7、命令。2 .同時繪制兩組數(shù)據(jù)的時序圖d=read.csv("double.csv",header=F)double=ts(d,start=1964,freq=1)plot(double,plot.type="multiple")#兩組數(shù)據(jù)兩個圖plot(double,plot.type="single")#兩組數(shù)據(jù)一個圖plot(double,plot.type="single",col=c("red","green"),lty=c(1,2)#設置每組數(shù)據(jù)圖的顏色、曲線類型)

8、3 .產(chǎn)生服從正態(tài)分布的隨機觀察值例題2.4隨機產(chǎn)生1000白噪聲序列觀察值d=rnorm(1000,0,1)#個數(shù)1000均值0方差1plot.ts(d)4 .純隨機性檢驗例題2.3續(xù)d=scan("temp.csv")temp=ts(d,freq=1,start=c(1949)Box.test(temp,type="Ljung-Box",lag=6)5 .差分計算x=1:10y=diff(x)k步差分k"xtxtk加入?yún)?shù)lag=k如計算x的3步差分為y=diff(x,lag=3)pp1p1p階差分xtxtxt1加入?yún)?shù)difference

9、s=p如2階差分僅XtXt1y=diff(x,differences=2)第三講例題3.1plot.ts(arima.sim(n=100,list(ar=0.8)# 模才AAR(1)模型，并作時序圖。plot.ts(arima.sim(n=100,list(ar=-1.1)# 非平穩(wěn)，無法得到時序圖。plot.ts(arima.sim(n=100,list(ar=c(1,-0.5)plot.ts(arima.sim(n=100,list(ar=c(1,0.5)例題3.5acf(arima.sim(n=100,list(ar=0.8)acf(arima.sim(n=100,list(ar=-1.

10、1)acf(arima.sim(n=100,list(ar=c(1,-0.5)acf(arima.sim(n=100,list(ar=c(1,0.5)例題3.7arima.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)例題2.5d=scan("a1.5.txt")#導入數(shù)據(jù)prop=ts(d,start=1950,freq=1)#轉化為時間序列數(shù)據(jù)plot(prop)#作時序圖a

11、cf(prop,12)#作自相關圖，拖尾pacf(prop,12)#作偏自相關圖，1階截尾Box.test(prop,type="Ljung-Box",lag=6)# 純隨機性檢驗,p值小于5%,序列為非白噪聲Box.test(prop,type="Ljung-Box",lag=12)arima(prop,order=c(1,0,0),method="ML")# 用AR(1膜型擬合，如參數(shù)method="CSS”,估計方法為條件最小二乘法，用條件最小二乘法時，不顯示AIC。arima(prop,order=c(1,0,0),

12、method="ML",include.mean=F)#用AR(1)模型擬合，不含截距項。tsdiag(arima(prop,order=c(1,0,0),method="ML")# 對估計進行診斷，判斷殘差是否為白噪聲summary(arima(prop,order=c(1,0,0),method="ML")a=arima(prop,order=c(1,0,0),method="ML")r=a$residuals#用r來保存殘差Box.test(r,type="Ljung-Box",lag=6

13、)#對殘差進行純隨機性檢驗predict(arima(prop,order=c(1,0,0),n.ahead=5)#預測未來5期prop.fore=predict(arima(prop,order=c(1,0,0),n.ahead=5)# 將未來5期預測值保存在prop.fore變量中U=prop.fore$pred+1.96*prop.fore$seL=prop.fore$pred-1.96*prop.fore$se#算出95%置信區(qū)間ts.plot(prop,prop.fore$pred,col=1:2)#作時序圖，含預測。lines(U,col="blue",lty=

14、"dashed")lines(L,col="blue",lty="dashed")#在時序圖中作出95%置信區(qū)間例題3.9d=scan("a1.22.txt")x=diff(d)arima(x,order=c(1,0,1),method="CSS")tsdiag(arima(x,order=c(1,0,1),method="CSS")第一點：對于第三講中的例2.5,運行命令arima(prop,order=c(1,0,0),method="ML")之后，顯

15、示：Call:arima(x=prop,order=c(1,0,0),method="ML")Coefficients:ar1intercept0.691481.5509s.e.0.09891.7453sigmaA2estimatedas15.51:loglikelihood=-137.02,aic=280.05注意：intercept下面的81.5509是均值，而不是截距！雖然intercept是截距的意思，這里如果用mean會更好。(themeanandtheinterceptarethesameonlywhenthereisnoARterm,均值和截距是相同的，只有在

16、沒有AR項的時候)如果想得到截距，利用公式計算。int=(1-0.6914)*81.5509=25.16661。課本P81的例2.5續(xù)中的截距25.17是正確的。第二點：如需計算參數(shù)的t統(tǒng)計量值和p值，利用下面的公式。ar的t統(tǒng)計量值=0.6914/0.0989=6.9909(注：數(shù)值與課本略有不同，因為課本用sas算的se=0.1029,R計算的se=0.0989)p值=pt(6.9909,df=48,lower.tail=F)*2pt()為求t分布求p值的函數(shù)，6.99為t統(tǒng)計量的絕對值，df為自由度=數(shù)據(jù)個數(shù)-參數(shù)個數(shù)，lower.tail=F表示所求p值為PT>t,如不加入這個參

17、數(shù)表示所求p值為PT<=t。乘2表示p值是雙側的(課本上的p值由sas算出，是雙側的)均值的t統(tǒng)計量值和p值同理。在時間序列中對參數(shù)顯著性的要求與回歸模型不同，我們更多的是考察模型整體的好壞，而不是參數(shù)。所以，R中的arima擬合結果中沒有給出參數(shù)的t統(tǒng)計量值和p值，如果題目沒有特別要求，一般不需要手動計算。第三點：修正第三講中的錯誤：例2.5中，我們用下面的語句對擬合arima模型之后的殘差進行了LB檢驗：a=arima(prop,order=c(1,0,0),method="ML")r=a$residualsa=arima(prop,order=c(1,0,0),

18、method="ML")r=a$residuals#用r來保存殘差Box.test(r,type="Ljung-Box",lag=6)#對殘差進行純隨機性檢驗最后一句不完整，需要加上參數(shù)fitdf=1,修改為Box.test(r,type="Ljung-Box",lag=6,fitdf=1)fitdf表示p+q,numberofdegreesoffreedomtobesubtractedifxisaseriesofresiduals,當檢驗的序列是殘差到時候，需要加上命令fitdf,表示減去的自由度。運行Box.test(r,type

19、="Ljung-Box",lag=6,fitdf=1)后，顯示的結果：Box.test(r,type="Ljung-Box",lag=6,fitdf=1)Box-Ljungtestdata:rX-squared=5.8661,df=5,p-value=0.3195“df=裱市自由度為5,由于參數(shù)lag=6,所以是滯后6期的檢驗。第四講#example4_1擬合線性模型x1=c(12.79,14.02,12.92,18.27,21.22,18.81,25.73,26.27,26.75,28.73,31.71,33.95)a=as.ts(x1)is.ts(a

20、)ts.plot(a)t=1:12tlm1=lm(at)summary(lm1)#返回擬合參數(shù)的統(tǒng)計量coef(lm1)#返回被估計的系數(shù)fitted(lm1)#返回模擬值residuals(lm1)#返回殘差值fit1=as.ts(fitted(lm1)ts.plot(a);lines(fit1,col="red")#擬合圖#eg1cs=ts(scan("eg1.txt",sep=",")csts.plot(cs)t=1:40lm2=lm(cst)summary(lm2)#返回擬合參數(shù)的統(tǒng)計量coef(lm2)#返回被估計的系數(shù)fi

21、t2=as.ts(fitted(lm2)#返回模擬值residuals(lm2)#返回殘差值ts.plot(cs);lines(fit2,col="red")#擬合圖#example4_2擬合非線性模型t=1:14x2=c(1.85,7.48,14.29,23.02,37.42,74.27,140.72,265.81,528.23,1040.27,2064.25,4113.73,8212.21,16405.95)x2plot(t,x2)m1=nls(x2-a*t+bAt,start=list(a=0.1,b=1.1),trace=T)summary(ml)#返回擬合參數(shù)的統(tǒng)

22、計量coef(m1)#返回被估計的系數(shù)fitted(ml)#返回模擬值residuals(ml)#返回殘差值plot(t,x2);lines(t,fitted(m1)#擬合圖#讀取excel中讀取文件，逗號分隔符a=read.csv("example4_2.csv",header=TRUE)t=a$tx=a$xxts.plot(x)m2=nls(xa*t+bAt,start=list(a=0.1,b=1.1),trace=T)summary(m2)#返回擬合參數(shù)的統(tǒng)計量coef(m2)#返回被估計的系數(shù)fitted(m2)#返回模擬值residuals(m2)#返回殘差值p

23、lot(t,x);lines(t,fitted(m2)#擬合圖#eg2I<-scan("eg2.txt")Ix=ts(data=I,start=c(1991,1),f=12)#化為時間序列xplot.ts(x)t=1:130t2=tA2m3=lm(xt+t2)coef(m3)#返回被估計的系數(shù)summary(m3)#返回擬合參數(shù)的統(tǒng)計量# 去不顯著的自變量，再次模擬m4=lm(xt2)coef(m4)#返回被估計的系數(shù)summary(m4)#返回擬合參數(shù)的統(tǒng)計量m2=fitted(m4)#返回模擬值y=ts(data=m2,start=c(1991,1),f=12)y

24、ts.plot(x);lines(y)# 平滑法# 簡單移動平均法x=c(5,5.4,5.8,6.2)xy=filter(x,rep(1/4,4),sides=1)y# 指數(shù)平滑for(iin1:3)x1=x1xi+1=0.25*xi+1+0.75*xi#HoltWintersFiltera=ts(read.csv("holt.csv",header=F),start=c(1978,1),f=1)am=HoltWinters(a,alpha=0.15,beta=0.1,gamma=FALSE,l.start=51259,b.start=4325)mfitted(m)plot

25、(m)plot(fitted(m)# 綜合cs=ts(read.csv("eg3.csv",header=F),start=c(1993,1),f=12)#讀取數(shù)據(jù)csts.plot(cs)#繪制時序圖cs.sea1=rep(0,12)cs.sea1for(iin1:12)for(jin1:8)cs.sea1i=cs.sea1i+csi+12*(j-1)cs.sea=(cs.sea/8)/(mean(cs)cs.seacs.sea2=rep(cs.sea,8)cs.sea2x=cs/cs.sea2xplot(x)t=1:96m1=lm(x-t)coef(m1)summary

26、(m1)m=ts(fitted(m1),start=c(1993,1),f=12)ts.plot(x,type="p");lines(m,col="red")r=residuals(m1)Box.test(r)#白噪聲檢驗第五講ffffffffffffffffffffffffffffffffffffffffffffffffTT，kk，kk，kk，TT，TTfjf>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>

27、;f>f>f># 回顧# 例5.1sha=ts(scan("sha.csv"),start=1964,freq=1)ts.plot(sha)diff(sha)par(mfrow=c(2,1)ts.plot(diff(sha)acf(diff(sha)# 例5.2car=ts(read.csv("car.csv",header=F),start=1950,freq=1)carpar(mfrow=c(3,1)ts.plot(car)ts.plot(diff(car)ts.plot(diff(car,differences=2)# 例5.3m

28、ilk=ts(scan("milk.txt"),start=c(1962,1),freq=12)milkpar(mfrow=c(3,1)ts.plot(milk)ts.plot(diff(milk)dm1=diff(diff(milk),lag=12)ts.plot(dm1)acf(dml)# 例5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)ffffffffffffffffffffffffffffffffffffffffffffffffffffff，kk，kk，f>f>f>f>f>f>f>f&g

29、t;f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f># 擬合ARIMA模型#5.8.1a=ts(scan("581.txt")par(mfrow=c(2,2)ts.plot(a)da=diff(a)ts.plot(da)acf(da,20)pacf(da,20)Box.test(da,6)fit1=arima(a,c(1,1,0),method="ML")predict(fit1,5)fffffffffffffff

30、fffffffffffffffffffffffffffffffffffffffffff，kk，f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>incom=ts(read.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(incom)ts.

31、plot(dincom)acf(dincom,lag=18)#自相關圖Box.test(dincom,type="Ljung-Box",lag=6)#白噪聲檢驗Box.test(dincom,type="Ljung-Box",lag=12)Box.test(dincom,type="Ljung-Box",lag=18)pacf(dincom,lag=18)fit1=arima(dincom,order=c(0,0,1),method="CSS")fit2=arima(incom,order=c(0,1,1),xre

32、g=1:length(incom),method="CSS")#見/stoffer/tsa2/Rissues.htmBox.test(fit2$resid,lag=6,type="Ljung-Box",fitdf=1)fore=predict(fit2,10,newxreg=(length(incom)+1):(length(incom)+10)# 疏系數(shù)模型# 例5.8w=ts(read.csv("w.csv"),start=1917,freq=1)w=w,1par(mfrow=c(2,

33、2)ts.plot(w)ts.plot(diff(w)acf(diff(w),lag=18)pacf(diff(w),lag=18)dw=diff(w)fit3=arima(dw,order=c(4,0,0),fixed=c(NA,0,0,NA,0),method="CSS")Box.test(fit3$resid,lag=6,type="Ljung-Box",fitdf=2)Box.test(fit3$resid,lag=12,type="Ljung-Box",fitdf=2)fit4=armaFit(arima(4,0,0),fi

34、xed=c(NA,0,0,NA),include.mean=F,data=dw,method="CSS")summary(fit4)# 例5.9ue=ts(scan("unemployment.txt"),start=1962,f=4)#讀取數(shù)據(jù)par(mfrow=c(2,2)#繪制時序圖ts.plot(ue)# 差分due=diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)# 平穩(wěn)性檢驗acf(ddue,lag=30)pacf(ddue,lag=30)arima(ddue,orde

35、r=c(0,0,0),method="ML")arima(ddue,order=c(4,0,0),method="ML")arma=arima(ddue,order=c(4,0,0),transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="ML")# 參數(shù)估計與檢驗(加載fArma程序包)fit2=armaFit(arima(4,0,0),include.mean=F,data=ddue,method="ML")summary(fit2)fit3=arm

36、aFit(arima(4,0,0),data=ddue,transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="CSS")summary(fit3)# 殘差白噪聲檢驗Box.test(arma$resid,6,fitdf=2,type="Ljung")# 擬合ft=ts(fitted(fit3),start=1963.25,f=4)dft=ts(rep(0,115),start=1963.25,f=4)for(iin1:115)dfti=fti+duei+uei+4ts.plot(ue);li

37、nes(dft,col="red")ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff“TTk，“TTk，TTk，kk，k，.，ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff# 例5.10乘積季節(jié)模型wue=ts(scan("wue.txt"),start=1948,f=12)arima(wue,order=c(1,1,1),seasonal=list(or

38、der=c(0,1,1),period=12),include.mean=F,method="CSS")ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff，kk，kk，”mrHIIIIIIffffffffffffffffffffffffffffffffffffffffffffffffffffffffffifff# 擬合Auto-Regressive模型eg1=ts(scan("582.txt")ts.plot(eg1)# 因變量關于時間的回歸模型fit.gls=gl

39、s(eg1-1+time(eg1),correlation=corARMA(p=1),method="ML")#seethenlmepackagesummary(fit.gls2)#theresults# 延遲因變量回歸模型leg1=lag(eg1,-1)y=cbind(eg1,leg1)fit=arima(y,1,c(0,0,0),xreg=y,2,include.mean=F)第六講# 回顧# 例5.1sha=ts(scan("sha.csv"),start=1964,freq=1)ts.plot(sha)diff(sha)par(mfrow=c(2

40、,1)ts.plot(diff(sha)acf(diff(sha)# 例5.2car=ts(read.csv("car.csv",header=F),start=1950,freq=1)carpar(mfrow=c(3,1)ts.plot(car)ts.plot(diff(car)ts.plot(diff(car,differences=2)# 例5.3milk=ts(scan("milk.txt"),start=c(1962,1),freq=12)milkpar(mfrow=c(3,1)ts.plot(milk)ts.plot(diff(milk)dm

41、1=diff(diff(milk),lag=12)ts.plot(dm1)acf(dm1)# 例5.5x=ts(cumsum(rnorm(1000,0,100)ts.plot(x)ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff# 擬合ARIMA模型# 上機指導5.8.1a=ts(scan("581.txt")par(mfrow=c(2,2)ts.plot(a)da=diff(a)ts.plot(da)acf

42、(da,20)pacf(da,20)Box.test(da,6)fit1=arima(a,c(1,1,0),method="ML")predict(fit1,5,newxreg=(length(a)+1):(length(a)+5)fit2=armaFit(arima(1,1,0),data=a,xreg=1:length(a),method="ML")summary(fit1)summary(fit2)#截距項不顯著，故舍去fit3=arima(a,c(1,1,0),method="ML")predict(fit3,5)ffffff

43、fffffffffffffffffffffffffffffffffffffffffffffffffffffjf>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f>f># 例5.8incom=ts(read.csv("incom.csv",header=F),start=1952,freq=1)incomts.plot(incom)dincom=diff(i

44、ncom)ts.plot(dincom)acf(dincom,lag=18)#自相關圖Box.test(dincom,type="Ljung-Box",lag=6)#白噪聲檢驗pacf(dincom,lag=18)fit=arima(incom,order=c(0,1,1),xreg=1:length(incom),method="CSS")# 見/stoffer/tsa2/Rissues.htmAutocorTest(fit$resid)#力口載FinTS包fore=predict(fit,10,newxr

45、eg=(length(incom)+1):(length(incom)+10)# 疏系數(shù)模型# 例5.8w=ts(read.csv("w.csv"),start=1917,freq=1)w=w,1par(mfrow=c(2,2)ts.plot(w)ts.plot(diff(w)acf(diff(w),lag=18)pacf(diff(w),lag=18)dw=diff(w)fit3=arima(dw,order=c(4,0,0),fixed=c(NA,0,0,NA,0),method="CSS")Box.test(fit3$resid,lag=6,typ

46、e="Ljung-Box",fitdf=2)Box.test(fit3$resid,lag=12,type="Ljung-Box",fitdf=2)fit4=armaFit(arima(4,0,0),fixed=c(NA,0,0,NA),include.mean=F,data=dw,method="CSS")#力口載fArma包，檢驗參數(shù)summary(fit4)# 例5.9# 讀取數(shù)據(jù)ue=ts(scan("unemployment.txt"),start=1962,f=4)# 繪制時序圖par(mfrow=c(

47、2,2)ts.plot(ue)# 差分due=diff(ue)ddue=diff(due,lag=4)ts.plot(ddue)Box.test(ddue,lag=6)# 平穩(wěn)性檢驗acf(ddue,lag=30)pacf(ddue,lag=30)arima(ddue,order=c(0,0,0),method="ML")arima(ddue,order=c(4,0,0),method="ML")arma=arima(ddue,order=c(4,0,0),transform.pars=F,fixed=c(NA,0,0,NA),include.mean=

48、F,method="ML")# 參數(shù)估計與檢驗(加載fArma程序包)fit2=armaFit(arima(4,0,0),include.mean=F,data=ddue,method="ML")summary(fit2)fit3=armaFit(arima(4,0,0),data=ddue,transform.pars=F,fixed=c(NA,0,0,NA),include.mean=F,method="CSS")summary(fit3)# 殘差白噪聲檢驗Box.test(arma$resid,6,fitdf=2,type="Ljung")# 擬合ft=ts(fitted(fit3),start=1963.25,f=4)dft=ts(rep(0,115),start=1963.25,f=4)for(iin1:115)dfti=fti+duei+