"Fossies" - the Fresh Open Source Software Archive

Member "ccmath-2.2.1/tseries/test/tfmest.c" (11 Dec 2000, 6235 Bytes) of package /linux/misc/old/ccmath-2.2.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 /*  tfmest.c    CCMATH mathematics library source code.
    2  *
    3  *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
    4  *  This code may be redistributed under the terms of the GNU library
    5  *  public license (LGPL). ( See the lgpl.license file for details.)
    6  * ------------------------------------------------------------------------
    7  */
    8 /*
    9     Test:  seqtsf  fixtsf
   10 
   11     Uses:  sdiff  evfmod  setevf  evmax  eigen
   12 
   13     Input parameters:  mod_fl -> name of file of model initialization data
   14 
   15                        ser_fl -> name of binary factor model time series
   16                                  data file [ created by gfarma ]
   17 
   18     Prompted input:  at prompt  ' (s,f) q->quit '
   19                           enter  s  for a sequential parameter update
   20                                  f  for a gauss-newton parameter update
   21                                  q  to terminate estimation
   22 
   23                      at prompt  ' save residuals? (y/n) '
   24                           enter  y  to save residuals to a file
   25                                  n  to exit without saving residuals
   26 */
   27 #include "ccmath.h"
   28 #include <math.h>
   29 struct mcof *pfc,*par,*pma;
   30 int nfc,nar,nma,np,ndif;
   31 struct fmod *x; int nmax;
   32 void main(int na,char **av)
   33 { double *var,*cr;
   34   struct mcof *pp;
   35   double ssq,sig,ev,*pe;
   36   int n,j,k; FILE *fm,*fs;
   37   char fnam[32],cs[4];
   38   if(na!=3){ printf("para: mod_fl ser_fl\n"); exit(-1);}
   39   fm=fopen(*++av,"r");
   40   printf("     model file: %s\n",*av);
   41 
   42 /* enter model initialization data */
   43   fscanf(fm,"%d %d %d %d",&nfc,&nar,&nma,&ndif);
   44   np=nfc+nar+nma;
   45 
   46 /* allocate store for madel parameters and variance */
   47   pfc=(struct mcof *)calloc(np,sizeof(*pfc));
   48   par=pfc+nfc; pma=par+nar;
   49   var=(double *)calloc(np*np+np,sizeof(*var));
   50   cr=var+np*np;
   51   for(j=0,pp=pfc; j<nfc ;++j,++pp){
   52     fscanf(fm,"%lf",&pp->cf);
   53    }
   54   for(; j<nfc+nar ;++j,++pp){
   55     fscanf(fm,"%d",&pp->lag); pp->lag-=1;
   56    }
   57   for(; j<np ;++j,++pp){
   58     fscanf(fm,"%d",&pp->lag); pp->lag-=1;
   59    }
   60   fclose(fm);
   61   printf(" Model Structure\n");
   62   printf(" %d initial factor parameters\n",nfc);
   63   for(j=0,pp=pfc; j<nfc ;++j) printf(" %.2f",(pp++)->cf);
   64   printf("\n");
   65   if(nar){
   66     printf(" %d ar para. at lags",nar);
   67     for(j=0,pp=par; j<nar ;++j) printf(" %d",(pp++)->lag+1);
   68     printf("\n");
   69    }
   70   if(nma){
   71     printf(" %d ma para. at lags",nma);
   72     for(j=0,pp=pma; j<nma ;++j) printf(" %d",(pp++)->lag+1);
   73     printf("\n");
   74    }
   75   printf("  difference order = %d\n",ndif);
   76 
   77 /* read factor model time series data */
   78   strcpy(fnam,*++av); fs=fopen(fnam,"rb");
   79   printf("     data file: %s\n",fnam);
   80   fread((void *)&nmax,sizeof(int),1,fs);
   81   x=(struct fmod *)calloc(nmax,sizeof(*x));
   82   n=fread((char *)x,sizeof(x[0]),nmax,fs);
   83 
   84 /* difference input series if required */
   85   if(ndif){
   86      for(j=k=0; j<n ;++j){
   87        ev=sdiff(x[j].val,ndif,j);
   88        if(j-ndif>=0){ x[k].val=ev; x[k++].fac=x[j].fac;}
   89       }
   90      n-=ndif;
   91    }
   92   printf(" %d points used in fit\n\n",n);
   93 
   94 /* start interactive estimation sequence */
   95   for(j=0; ;++j){
   96     fprintf(stderr," (s/f) q->quit ");
   97     if(*gets(cs)=='q') break;
   98 
   99 /* sequential estimation step */
  100     if(cs[0]=='s') ssq=seqtsf(x,n,var,j);
  101 
  102 /* gauss-newton estimation step */
  103     else ssq=fixtsf(x,n,var,cr);
  104 
  105     printf("%d  ssq= %8.3f  ",j,ssq);
  106     if(cs[0]=='s') printf("seq.\n"); else printf("fix\n");
  107     printf(" p_vec: ");
  108     for(pp=pfc,k=0; k<np ;++k,++pp) printf("  %7.4f",pp->cf);
  109 /* compute maximum eigenvalue and vector of parameter variance */
  110     ev=evmax(var,cr,np)*sqrt(ssq/(n-np));
  111     printf("\n max ev and vector : %9.6f\n",ev);
  112     matprt(cr,1,np," %.3f");
  113     printf("\n");
  114    }
  115   printf(" final state:\n");
  116   printf("  ssq = %8.3f\n",ssq);
  117   printf("  p_vec: ");
  118   for(pp=pfc,k=0; k<np ;++k,++pp) printf("  %7.4f",pp->cf);
  119   printf("\n");  
  120   sig=sqrt(ssq/(n-np)); printf("  rms e = %8.5f\n",sig);
  121   for(k=0; k<np*np ;) var[k++]*=sig;
  122   printf("  variance matrix:\n"); matprt(var,np,np," %10.7f");
  123 
  124 /* compute all eigenvalues and vectors of final variance estimate */
  125   eigen(var,cr,np);
  126   printf("  eigenvalues:\n"); matprt(cr,1,np," %10.7f");
  127   printf("  vectors:\n"); matprt(var,np,np," %10.3f");
  128 
  129 /* compute and save the model residuals in a file with .er extension */
  130   fprintf(stderr," save model residuals? (y/n) "); gets(cs);
  131   if(cs[0]=='y'){
  132     for(j=0; fnam[j] && fnam[j]!='.' ;++j);
  133     fnam[j]=0; strcat(fnam,".er");
  134     fm=fopen(fnam,"wb");
  135     pe=(double *)calloc(n,sizeof(*pe));
  136     setevf(1);
  137     for(j=0; j<n ;++j) pe[j]=evfmod(x[j]);
  138 /* write header specifing number of residuals to file */
  139     fwrite((void *)&n,sizeof(int),1,fm);
  140 /* write residuals */
  141     fwrite((char *)pe,sizeof(double),n,fm);
  142    }
  143 }
  144 /*  Test output
  145 
  146      model file: data/tfs0.d
  147  Model Structure
  148  2 initial factor parameters
  149  0.48 0.48
  150  2 ar para. at lags 1 2
  151  1 ma para. at lags 1
  152   difference order = 0
  153      data file: data/tfs0.b
  154  400 points used in fit
  155 
  156 0  ssq=  404.842  seq.
  157  p_vec:   -0.2561   0.7800   0.8560  -0.4883  -0.4139
  158  max ev and vector :  0.022466
  159  0.706 0.708 0.001 -0.010 0.008
  160 
  161 1  ssq=  367.338  seq.
  162  p_vec:   -0.2520   0.7859   0.8541  -0.4915  -0.4384
  163  max ev and vector :  0.011396
  164  0.707 0.708 0.001 -0.007 0.006
  165 
  166 2  ssq=  365.777  seq.
  167  p_vec:   -0.2504   0.7879   0.8531  -0.4924  -0.4483
  168  max ev and vector :  0.007779
  169  0.707 0.708 0.000 -0.005 0.005
  170 
  171 3  ssq=  364.058  fix
  172  p_vec:   -0.2463   0.7912   0.8472  -0.4919  -0.4793
  173  max ev and vector :  0.024584
  174  0.707 0.708 -0.000 -0.001 0.001
  175 
  176  final state:
  177   ssq =  364.058
  178   p_vec:   -0.2463   0.7912   0.8472  -0.4919  -0.4793
  179   rms e =  0.96003
  180   variance matrix:
  181   0.0128710  0.0116984  0.0000696 -0.0000756  0.0001102
  182   0.0116984  0.0128999 -0.0000844  0.0000451 -0.0000749
  183   0.0000696 -0.0000844  0.0035962 -0.0024481  0.0026019
  184  -0.0000756  0.0000451 -0.0024481  0.0030263 -0.0021005
  185   0.0001102 -0.0000749  0.0026019 -0.0021005  0.0041369
  186   eigenvalues:
  187   0.0245840  0.0084018  0.0011812  0.0015415  0.0008221
  188   vectors:
  189       0.707      0.018      0.706     -0.036      0.005
  190       0.708     -0.019     -0.705      0.037     -0.006
  191      -0.000      0.596      0.008      0.336     -0.729
  192      -0.001     -0.513     -0.009     -0.539     -0.668
  193       0.001      0.617     -0.057     -0.770      0.149
  194 */