nonlin_singleprop2.r

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

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

M<-200 ##number of particles
storage.mode(M) <- "integer"
storage.mode(n) <- "integer"
storage.mode(z)<-"double"
storage.mode(ss)<-"double"
storage.mode(sp)<-"double"
storage.mode(sn)<-"double"
UP<-as.integer(10000/M)

DF <- 5 ##degrees of freedom for T-proposal
storage.mode(DF) <- "double"
 

N<-100 ##number of independen analyses
##update progress every UP iterations
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)
ests.st<-matrix(0,nrow=N,ncol=n)
pr.st<-matrix(0,nrow=N,ncol=n)
LEst<-rep(0,N)
for(j in 1:N){

  LE<-0 #log-evidence
  est<-rep(0,n) #Filtered Estimates
  ests<-rep(0,n) #Smoothed Estimates
  pr.est <- rep(0,n) ##Filter estimates of Pr(alpha<true_alpha)
  
  storage.mode(LE)<-"double"
  storage.mode(est)<-"double"
  storage.mode(ests)<-"double"
  storage.mode(pr.est)<-"double"

  val<-rep(0,n)
  storage.mode(val) <- "double"
  SMOOTH<-0 ##do smoothing 0=No; 1=Yes
  OUTPUT<-1 ##output particles and weights?
  storage.mode(SMOOTH) <- "integer"
  storage.mode(OUTPUT) <- "integer"

  out <- .C("NL_singleprop2",z,n,M,sn,ss,sp,DF,SMOOTH,val,LE=LE,est=est,ests=ests,pr.est=pr.est,OUTPUT)

 

if(as.integer(j/UP)==j/UP) cat("End of It ",j,"\n")
LEst[j]<-out$LE
est.st[j,]<-out$est
ests.st[j,]<-out$ests
pr.st[j,] <- out$pr.est
}

##Print Results
cat("Var ",var(LEst),"\n")
cat("Mean ",mean(LEst),"\n")

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


####If output ==1 then read in weighted particles and plot density estiamtes

#scan

part.st<-matrix(scan("NLsp_PARTICLES.txt"),nrow=n,ncol=M,byrow=T)
w.st<-matrix(scan("NLsp_WEIGHTS.txt"),nrow=n,ncol=M,byrow=T)

#plot
par(mfrow=c(4,3))
for(t in 1:10){
  index<-sample(1:M,size=M,prob=w.st[t,],replace=T);  part<-part.st[t,index];  plot(density(part),type="l"); lines(c(alpha[t],alpha[t]),c(0,1),col=2);
}