解禁します

計量時系列分析の課題提出日を過ぎたので、書いたコードを解禁します。

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)