arch_singleprop.r

用MATLAC做的非线性时间序列预测和平滑,可以用来做非线性模型的预测和估计

#Load in .so file (C code)
dyn.load("ARCH_singleprop.so")


#number of particles K resampled from M
M<-200
K<-200
storage.mode(M)<-"integer"
storage.mode(K)<-"integer"
storage.mode(n)<-"integer"
storage.mode(z)<-"double"
storage.mode(sig)<-"double"
storage.mode(b)<-"double"

UP<-as.integer(10000/M)
             
               
RES<-0 #IGNORE
SMOOTH <- 1 #SMOOTH (0=No; 1=Yes. 
df<-4 #df in t-proposal
FULL <- K #DO NOT CHANGE!
##RES AND FULL WHERE INITIALLY INTRODUCED TO ALLOW FOR DIFFERENT VERSIONS FOR
##CALCULATING THE WEIGHTS/RESAMPLING. USING FULL=K PRODUCE RPF IN PAPER - FOR
##THIS RES-VALUE SHOULD NOT MATTER.
storage.mode(RES)<-"integer"
storage.mode(FULL)<-"integer"
storage.mode(df)<-"integer"


#number of simulations
N<-100

###This is just to print out progress
if(UP>as.integer(N/10)) UP<-as.integer(N/10)
if(UP==0) UP <- 1

###storage of results
est.st<-matrix(0,nrow=N,ncol=n)
pr.st<-matrix(0,nrow=N,ncol=n)
est2.st<-matrix(0,nrow=N,ncol=n)
LEst<-rep(0,N)


for(j in 1:N){
  ###This just samples QMC from the prior distribution
  u<-runif(K,0,1/K)+(0:(K-1))/K 
  init <- qnorm(u,0,1)
  storage.mode(init) <- "double"
  LE<-0 #log-evidence (up to some additive constant!)
  est<-rep(0,n) #Filter estimate of theta_1,...,theta_n
  pr.est <- rep(0,n) #Filter estimate of pr(theta_t< theta_true_t)
  est2<-rep(0,n) #Smooth estimate of theta_1,...,theta_n
  storage.mode(LE)<-"double"
  storage.mode(est)<-"double"
  storage.mode(pr.est)<-"double"
  storage.mode(est2)<-"double"
  storage.mode(theta) <- "double"
  ECHO <- (j==1) ##Give details of progress?
  storage.mode(ECHO)<-"integer"
  storage.mode(SMOOTH)<-"integer"

  out <- .C("ARCH_singleprop",z,n,M,K,sig,b,FULL,RES,df,ECHO,SMOOTH,init,theta,LE=LE,est=est,est2=est2,pr.est=pr.est)


if(as.integer(j/UP)==j/UP) cat("End of It ",j,"\n")
LEst[j]<-out$LE
est.st[j,]<-out$est
est2.st[j,]<-out$est2
pr.st[j,] <- out$pr.est
}
####Print results of Simulations
cat("Var ",var(LEst),"\n")
cat("Mean ",mean(LEst),"\n")

cat("Est ", mean(apply(est.st,2,var)),"\n")
cat("Est2 ",mean(apply(est2.st,2,var)),"\n")
cat("Pr ",mean(apply(pr.st,2,var)),"\n")