計量時系列分析の課題提出日を過ぎたので、書いたコードを解禁します。
ts.work<-function(gamma0,phi,n=500,no=c("kadai1","kadai2","kadai3","kadai4")){ no<-match.arg(no) ts.work<-list(gamma0=gamma0,phi=phi,n=n,no=no) class(ts.work)<-"ts.work" return(ts.work) } gamma.ts.work<-function(ts.work){ gamma0<-ts.work$gamma0 phi<-ts.work$phi n<-ts.work$n no<-ts.work$no do<-no x<-c() kadai1<-function(){ x[1]<-rnorm(1)*sqrt(gamma0) for(i in 2:n){ x[i]<-phi[1]*x[i-1]+rnorm(1) } return(x) } kadai2and3<-function(){ x[1]<-rnorm(1)*sqrt(gamma0) x[2]<-rnorm(1)*sqrt(gamma0) for(i in 3:n){ x[i]<-phi[1]*x[i-1]+phi[2]*x[i-2]+rnorm(1) } return(x) } kadai4<-function(){ x[1]<-rnorm(1)*sqrt(gamma0) x[2]<-rnorm(1)*sqrt(gamma0) x[3]<-rnorm(1)*sqrt(gamma0) for(i in 4:n){ x[i]<-phi[1]*x[i-1]+phi[2]*x[i-2]+phi[3]*x[i-3]+rnorm(1) } return(x) } switch(do, kadai1 = kadai1(), kadai2 = kadai2and3(), kadai3 = kadai2and3(), kadai4 = kadai4() ) } rho<-function(ts.work){ rho<-function(ts.work,n){ gamma0<-ts.work$gamma0 phi<-ts.work$phi no<-ts.work$no do<-no n<-n rho1<-function(n){ return(0.75^(n)) } rho2<-function(n) { if (n == 0) return(1) else if(n==1) return(0.25+0.125*0.25/(1-0.125)) else return(0.25*Recall(n-1)+0.125*Recall(n-2)) } rho3<-function(n) { if (n == 0) return(1) else if(n==1) return(1-0.5*1/(1+0.5)) else return(Recall(n-1)-0.5*Recall(n-2)) } rho4<-function(n){ if (n == 0) return(1) else if(n==1) return((phi[1]+phi[2]*phi[3])/(1-phi[2]-phi[1]*phi[3]-phi[3]^2)) else if(n==2) return(((phi[1]+phi[3])*(phi[1]+phi[2]*phi[3])) /(1-phi[2]-phi[1]*phi[3]-phi[3]^2)+phi[2]) else return(1.5*Recall(n-1)-1*Recall(n-2)+0.25*Recall(n-3)) } switch(do, kadai1 = rho1(n), kadai2 = rho2(n), kadai3 = rho3(n), kadai4 = rho4(n)) } x<-c() for(i in 1:21){ x[i]<-rho(ts.work,i-1) } return(x) } plot.ts.work<-function(ts.work){ x<-gamma(ts.work) no<-ts.work$no png(paste("kadai",substring(no,6,6),"_ts1.jpg",sep="")) |p!| plot(x) title("元のデータの時系列プロット") dev.off() png(paste("kadai",substring(no,6,6),"_ts2.jpg",sep="")) plot(x) lines(x) title("元のデータの時系列プロット") dev.off() theory.plot<-function(ts.work,n){ x<-rho(ts.work) png(paste("kadai",substring(no,6,6),"_theory",n,".jpg",sep="")) plot(0:20,acf(gamma(ts.work),plot=FALSE)$acf[1:21,,] ,xlim=c(0,20),ylim=c(-0.2,1),xlab="時間間隔",ylab="自己相関係数") abline(0,0) title("標本自己相関係数と理論自己相関係数の比較") lines(0:20,x) dev.off() } for(i in 1:5){ theory.plot(ts.work,i) } } ts <- ts.work(gamma0=1/(1-9/16),phi=0.75,no="kadai1") gamma(ts) acf(gamma(ts),plot=FALSE) rho(ts) plot(ts) ts <- ts.work(gamma0=448/405,phi=c(0.25,0.125),no="kadai2") gamma(ts) acf(gamma(ts),plot=FALSE) rho(ts) plot(ts) ts <- ts.work(gamma0=12/5,phi=c(1,-0.5),no="kadai3") gamma(ts) acf(gamma(ts),plot=FALSE) rho(ts) plot(ts) ts <- ts.work(gamma0=as.vector(solve(matrix(c(4,-7,4,6,-8,1,37,-56,12),3,3,byrow=TRUE)) %*%c(0,0,-16))[1],phi=c(1.5,-1,0.25),no="kadai4") gamma(ts) acf(gamma(ts),plot=FALSE) rho(ts) plot(ts)