arch_singleprop.c

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

#include <math.h>
#include "rand.c"

/*analyse n observations from the ARCH model
parameters are sig and b; K particles subsampled from M

M order particles are obtained via generating from a single proposal- t density with df degrees of freedom

output log-evidence

Full determines weight calculated 
0 = based on sample of 1 current particle
>0 = based on average over all current particles

Res determines resampling algorithm



Prior for theta[1] is rnorm(0,1); init are stratified sample from N(0,1)
*/
double qgen(double *qt,long Nq, double U);
void FCinit(double *w,int M,int O, double *c,int *n1,int *n2,int *in1,int *in2);
double Pr(double x2, double x1, double *b);
void ARCH_singleprop(double *z,long *n,long *M,long *K,double *sig,double *b,long *FULL,long *RES, long *df, long *ECHO, long *SMOOTH, double *init, double *val, double *LE, double *est,double *est2,double *pr_est)
{

  double *weight,*weight_new,*Ust;
  double **part_st,**w_st;
  double *w_temp,*w2;
  int *in1,*in2;
  double c;
  int n1,n2;
  int *index,temp;
  double *particles, *v,*w,*mu,*sd,*particles_new;
  long i,j,t;
  double sum,sig2,U;
  double Me,V;
  double xlast,plast,xnew,pnew,qinv;
  double *qt;
  long Nq;

  FILE *seed_file;
 
  /*obtain seed*/
  seed_file = fopen(".seed","r");
  if (seed_file==NULL){
    (void) fprintf(stderr,"Can not open .seed\n Using default values.");
    *seed=1349761;
    *(seed+1)=2498265;
  }else  {
    seed_set(seed_file);
    (void)fclose(seed_file);
  }
            
  /*look-up table*/
  Nq=5* *M;
  if(Nq<1000) Nq=1000;
  qt=(double *)calloc(Nq,sizeof(double));

  /*calculate look-up table*/
  i=Nq-1;
  plast=0.5;
  xlast=0.0;
  pnew=0.5;
  xnew=0.0;
  while(i>=0){
    if(pnew< (i+0.5)/(2.0*Nq)){
      if(pnew==plast)
	qt[i]=xnew;
      else{
	qt[i]=xnew+ ((i+0.5)/(2.0*Nq)-pnew)*(xlast-xnew)/(plast-pnew);
      }
      i--;
    }else{
      plast=pnew;
      xlast=xnew;
      xnew=xlast-0.001;
      pnew=pt((double) *df,xnew);
    }
  }
  
  if(*ECHO)
  (void)fprintf(stderr,"Finished calculating Look-up table\n");


  sig2=*sig * *sig;

  if(*FULL>0){
    in1=(int *)calloc(*FULL,sizeof(int));
  }
  in2=(int *)calloc(*M,sizeof(int));

  part_st=(double **)calloc(*n,sizeof(double *));
  w_st=(double **)calloc(*n,sizeof(double *));
  for(i=0;i< *n;i++){
    part_st[i]=(double *)calloc(*K,sizeof(double));
    w_st[i]=(double *)calloc(*K,sizeof(double));
  }

  particles=(double *)calloc(*K,sizeof(double));
  weight=(double *)calloc(*K,sizeof(double));
  weight_new=(double *)calloc(*M,sizeof(double));
  Ust=(double *)calloc(*M,sizeof(double));
  index = (int *)calloc(*M,sizeof(int));
  particles_new=(double *)calloc(*M,sizeof(double));
  v=(double *)calloc(*K,sizeof(double));
  w=(double *)calloc(*K,sizeof(double));
  w_temp=(double *)calloc(*K,sizeof(double));
  w2=(double *)calloc(*K,sizeof(double));
  mu=(double *)calloc(*K,sizeof(double));
  sd=(double *)calloc(*K,sizeof(double));

  /*initiate*/
  mu[0]=z[0]/(1+ sig2);
  sd[0]= sqrt(sig2/(1+sig2));
  for(i=0;i<*K;i++){
    particles[i]=mu[0]+sd[0]*init[i];
    part_st[0][i]=particles[i];
    weight[i]=1/(double)*K;
    w_st[0][i]=weight[i];
  }
  *LE= -0.5*z[0]*z[0]/(sig2+1.0) -0.5*log(sig2 +1.0);
  est[0]=0.0;
  pr_est[0]=0.0;
  for(i=0;i<*K;i++){
    est[0]+= particles[i]*weight[i];
    if(particles[i]<val[0])
      pr_est[0]+=weight[i];
  }
  for(t=1;t<*n;t++){
    sum=0.0;
    /*calculate proposal for time t*/
    for(i=0;i<*K;i++){
      v[i]=b[0]+b[1]*particles[i]*particles[i];
      w[i]=exp(-0.5*(z[t]*z[t])/(sig2 +v[i]) -0.5*log(v[i]+sig2))*weight[i];
      sum+=w[i];
      mu[i]=v[i]*z[t]/(sig2+v[i]);
      sd[i]=sqrt(v[i]/(v[i]+sig2))* *sig;
    }
    *LE += log(sum);
    /*calculate approximation - proposal at next time step*/
    Me=0.0;
    V=0.0;
    for(i=0;i< *K;i++)
      Me+=w[i]*mu[i];
    Me/=sum;
    for(i=0;i< *K;i++){
      V+= w[i]*(sd[i]*sd[i]+(mu[i]-Me)*(mu[i]-Me));
    }
    V/=sum;
    V=exp(0.5*log(V));
    /*sample *M new particles and weight*/
    /*U=runif();
      U/= *M;*/
    for(i=0;i< *M;i++){
      U= (runif()+(double)i)/(double)*M;
      particles_new[i]= qgen(qt,Nq,U);
      particles_new[i]=particles_new[i]*V+Me;
      /* U+= 1/(double)*M;*/
    }
    if(*FULL>0){
      if(*FULL>= *M){
	for(i=0;i<*M;i++){
	  weight_new[i]=0.0;
	  qinv=1.0/dt((double) *df,(particles_new[i]-Me)/V);
	  for(j=0;j<*K;j++)
	    weight_new[i]+= w[j]*dnorm((particles_new[i]-mu[j])/sd[j])/sd[j];
	  weight_new[i]*=qinv;
	}
      }else{
	/*calculate c,n1,n2,in1,in2 etc*/
	FCinit(w,*M,*FULL,&c,&n1,&n2,in1,in2);
	for(i=0;i<*M;i++){
	  weight_new[i]=0.0;
	  qinv=1.0/dt((double) *df,(particles_new[i]-Me)/V);
	  for(j=0;j<n1;j++)
	    weight_new[i] += w[in1[j]]*dnorm((particles_new[i]-mu[in1[j]])/sd[in1[j]]); 
	  if(n1< *FULL){
	    /*sample from in2*/
	    U=sum*runif()/ ((double)*FULL-n1);
	    j=0;
	    while(j < n2){
	      U-= w[in2[j]];
	      if(U<=0){
		weight_new[i]+= c* dnorm((particles_new[i]-mu[in2[j]])/sd[in2[j]]); 
		U+= sum/((double)*FULL-n1);
	      }
	      j++;
	    }

	  }
	  weight_new[i]*=qinv; 
	  	
	}
      }

    }else{
      /*sample index*/
      for(i=0;i<*M;i++)
	Ust[i]= (runif()+(double)i)/(double)*M;
      rdunif(w,sum,index,Ust,*M,*K);
       /*permute index*/
      for(i=0;i<*M;i++){
	j=(int)(runif()* (*M-i));
	temp=index[j];
	index[j]=index[*M-i-1];
	index[*M-i-1]=temp;
      }
      for(i=0;i<*M;i++){
	weight_new[i]= dnorm((particles_new[i]-mu[index[i]])/sd[index[i]]); /*true prob -weights cancel*/
	weight_new[i]/= dt((double) *df,(particles_new[i]-Me)/V); /*proposal prob*/
      }
    }
    /*normalise weights*/
    sum=0.0;
    for(i=0;i< *M;i++){
      sum+=weight_new[i];
    }
    /**LE+= log(sum/ *M);*/
    for(i=0;i< *M;i++) weight_new[i]=weight_new[i]/sum;

    /*resample new_particles*/
    /*sample *M new particles*/
    if(*K!=*M){
      if(*RES<=0){ /*simple resampling*/
	j=0;
	i=0;
	U=runif();
	U*= 1/(double)*K;
	while(i<*M){
	  if(U<= 1/(double)*M){
	    particles[j]=particles_new[i];
	    j++;
	    U+= 1/(double)*K;
	  }else{
	    i++;
	    U-=1/(double)*M;
	  }
	}
	if(j!= *K)
	  (void)fprintf(stderr,"Error - wrong number of particles sampled \n");
      }else{
	/*FC resampling*/
	
      }
    }else{
      for(i=0;i<*K;i++){
	particles[i]=particles_new[i];
	weight[i]=weight_new[i];
      }
    }
    est[t]=0.0;
    pr_est[t]=0.0;
    for(i=0;i<*K;i++){
      if(particles[i]<val[t])
	pr_est[t]+=weight[i];
      est[t]+= particles[i]*weight[i];
      w_st[t][i]=weight[i];
      part_st[t][i]=particles[i];
    }


  }

  if(*SMOOTH){
    /*smoothing - marginal*/
    est2[*n-1]=est[*n-1];
    for(t= *n-2;t>=0;t--){
      for(i=0;i< *K;i++)
	w2[i]=0.0;

      for(j=0;j<*K;j++){
	sum=0.0;
	for(i=0;i< *K;i++){
	  w_temp[i]=w_st[t][i]*Pr(part_st[t+1][j],part_st[t][i],b);
	  sum+=w_temp[i];
	}
	for(i=0;i< *K;i++)
	  w2[i]+= w_temp[i]*weight[j]/sum;
      }
      est2[t]=0.0;
      for(i=0;i< *K;i++){
	est2[t]+= w2[i] * part_st[t][i];
	weight[i]=w2[i];
      }
    }
  }


  /*change .seed file*/
  seed_file=fopen(".seed","w");
  seed_put(seed_file);
  (void)fclose(seed_file);
  

  free(Ust);
  free(index);
  free(qt);
  free(sd);
  free(particles);
  free(weight);
  free(weight_new);
  free(particles_new);
  free(mu);
  free(v);
  free(w);
  free(w2);
  free(w_temp);
  for(i=0;i < *n;i++){
    free(w_st[i]);
    free(part_st[i]);
  }
  
  free(w_st);
  free(part_st);
 

  if(*FULL>0){
    free(in1);
  }
  free(in2);



}

/*
calculate the Uth quantile given

qt[i]=x then Pr(t<x)= (1+0.5)/(2*Nq)
*/
double qgen(double *qt,long Nq, double U)
{
  double p;
  long i;
  if(U<=0.5){
    if(U> 0.5- 0.25/(double)Nq){
      return( qt[Nq-1]*(0.5-U)*4*Nq );
    }
    if(U< 0.25/(double)(Nq)){
      p=0.25/(double)(Nq);
      return( qt[0]+(U-p)*(qt[1]-qt[0])/(2.0*p) );
    }
    i=(int)(2*Nq*U+0.5);
    p= (i-0.5)/(2.0*Nq);
    return( qt[i-1]+(U-p)*(qt[i]-qt[i-1])*2.0*Nq);

  }else{
    U=1-U;
    if(U> 0.5- 0.25/(double)Nq){
      return( -qt[Nq-1]*(0.5-U)*4*Nq );
    }
    if(U< 0.25/(double)(Nq)){
      p=0.25/(double)(Nq);
      return( -qt[0]-(U-p)*(qt[1]-qt[0])/(2.0*p) );
    }
    i=(int)(2*Nq*U+0.5);
    p= (i-0.5)/(2.0*Nq);
    return( -qt[i-1]-(U-p)*(qt[i]-qt[i-1])*2.0*Nq);
 

  }

}
/*initiate for the FC sampling
input is w: weights
M: number of weights
N: number of samples

output
in1 (length n1): weights to keep;
in2 (length n2): others
c : resampled weight
*/
void FCinit(double *w,int M,int N, double *c,int *n1,int *n2,int *in1,int *in2)
{
  int i,j;
  double sw,w1;
  int n1old;


  /* calculate sum w*/
  sw=0.0;
  for(i=0;i<M;i++)
    sw+=w[i];
  /*initiate*/
  *n1=0;
  n1old= -1;
  *n2=M;
  for(i=0;i<M;i++)
    in2[i]=i;
  w1=0;
  *c= sw/ (double)N;

  while(*n1 != n1old && w1<sw && *n1<N){
    n1old= *n1;
    for(i=0;i< *n2;i++){
      if(w[in2[i]] >= *c){
	in1[*n1]=in2[i];
	*n1+=1;
	for(j=i+1;j< *n2;j++)
	  in2[j-1]=in2[j];
	*n2-=1;
	w1+=w[in2[i]];
      }
    }
    if(*n1<N)
      *c = (sw-w1)/((double)(N-*n1));
    else
      *c=0.0;
  }


}

/*prob of moving x1 to x2*/
double Pr(double x2, double x1, double *b)
{
  double z;
  double var,sd;

  var=b[0]+b[1]*x1*x1;
  sd=sqrt(var);
  z=x2/sd;
  return(dnorm(z)/sd );


}