dfa.c

消除趋势分析（DFA）目前已经成为时间序列分析的重要方法

    long i;
    int j, k;

    beta = vector(1, nfit);
    covar = matrix(1, nfit, 1, nfit);
    covar0 = matrix(1, nfit, 1, nfit);
    indxc = ivector(1, nfit);
    indxr = ivector(1, nfit);
    ipiv = ivector(1, nfit);
    mse = vector(1, nr);
    x = matrix(1, rs[nr], 1, nfit);
    for (i = 1; i <= rs[nr]; i++) {
	x[i][1] = 1.0;
	x[i][2] = i;
	for (j = 3; j <= nfit; j++)
	    x[i][j] = x[i][j-1] * i;
    }
}

/* This function frees all memory previously allocated by this program. */
void cleanup()
{
    free_matrix(x, 1, rs[nr], 1, nfit);
    free_vector(mse, 1, nr);
    free_ivector(ipiv, 1, nfit);
    free_ivector(indxr, 1, nfit);
    free_ivector(indxc, 1, nfit);
    free_matrix(covar0, 1, nfit, 1, nfit);
    free_matrix(covar, 1, nfit, 1, nfit);
    free_vector(beta, 1, nfit);
    free_lvector(rs, 1, rslen);	/* allocated by rscale() */
    free(seq);			/* allocated by input() */
}

static char *help_strings[] = {
 "usage: %s [OPTIONS ...]\n",
 "where OPTIONS may include:",
 " -d K             detrend using a polynomial of degree K",
 "                   (default: K=1 -- linear detrending)",
 " -h               print this usage summary",
 " -i               input series is already integrated",
 " -l MINBOX        smallest box width (default: 2K+2)",
 " -s               sliding window DFA",
 " -u MAXBOX        largest box width (default: NPTS/4)",
 "The standard input should contain one column of data in text format.",
 "The standard output is two columns: log(n) and log(F) [base 10 logarithms],",
 "where n is the box size and F is the root mean square fluctuation.",
NULL
};

void help(void)
{
    int i;

    (void)fprintf(stderr, help_strings[0], pname);
    for (i = 1; help_strings[i] != NULL; i++)
	(void)fprintf(stderr, "%s\n", help_strings[i]);
}

/* polyfit() is based on lfit() and gaussj() from Numerical Recipes in C
   (Press, Teukolsky, Vetterling, and Flannery; Cambridge U. Press, 1992).  It
   fits a polynomial of degree (nfit-1) to a set of boxsize points given by
   x[1...boxsize][2] and y[1...boxsize].  The return value is the sum of the
   squared errors (chisq) between the (x,y) pairs and the fitted polynomial.
*/
double polyfit(double **x, double *y, long boxsize, int nfit)
{
    int icol, irow, j, k;
    double big, chisq, pivinv, temp;
    long i;
    static long pboxsize = 0L;

    /* This block sets up the covariance matrix.  Provided that boxsize
       never decreases (which is true in this case), covar0 can be calculated
       incrementally from the previous value. */
    if (pboxsize != boxsize) {	/* this will be false most of the time */
	if (pboxsize > boxsize)	/* this should never happen */
	    pboxsize = 0L;
	if (pboxsize == 0L)	/* this should be true the first time only */
	    for (j = 1; j <= nfit; j++)
		for (k = 1; k <= nfit; k++)
		    covar0[j][k] = 0.0;
	for (i = pboxsize+1; i <= boxsize; i++)
	    for (j = 1; j <= nfit; j++)
		for (k = 1, temp = x[i][j]; k <= j; k++)
		    covar0[j][k] += temp * x[i][k];
	for (j = 2; j <= nfit; j++)
	    for (k = 1; k < j; k++)
		covar0[k][j] = covar0[j][k];
	pboxsize = boxsize;
    }
    for (j = 1; j <= nfit; j++) {
	beta[j] = ipiv[j] = 0;
	for (k = 1; k <= nfit; k++)
	    covar[j][k] = covar0[j][k];
    }
    for (i = 1; i <= boxsize; i++) {
	beta[1] += (temp = y[i]);
	beta[2] += temp * i;
    }
    if (nfit > 2)
	for (i = 1; i <= boxsize; i++)
	    for (j = 3, temp = y[i]; j <= nfit; j++)
		beta[j] += temp * x[i][j];
    for (i = 1; i <= nfit; i++) {
	big = 0.0;
	for (j = 1; j <= nfit; j++)
	    if (ipiv[j] != 1)
		for (k = 1; k <= nfit; k++) {
		    if (ipiv[k] == 0) {
			if ((temp = covar[j][k]) >= big ||
			    (temp = -temp) >= big) {
			    big = temp;
			    irow = j;
			    icol = k;
			}
		    }
		    else if (ipiv[k] > 1)
			error("singular matrix");
		}
	++(ipiv[icol]);
	if (irow != icol) {
	    for (j = 1; j <= nfit; j++) SWAP(covar[irow][j], covar[icol][j]);
	    SWAP(beta[irow], beta[icol]);
	}
	indxr[i] = irow;
	indxc[i] = icol;
	if (covar[icol][icol] == 0.0) error("singular matrix");
	pivinv = 1.0 / covar[icol][icol];
	covar[icol][icol] = 1.0;
	for (j = 1; j <= nfit; j++) covar[icol][j] *= pivinv;
	beta[icol] *= pivinv;
	for (j = 1; j <= nfit; j++)
	    if (j != icol) {
		temp = covar[j][icol];
		covar[j][icol] = 0.0;
		for (k = 1; k <= nfit; k++) covar[j][k] -= covar[icol][k]*temp;
		beta[j] -= beta[icol] * temp;
	    }
    }
    chisq = 0.0;
    if (nfit <= 2)
	for (i = 1; i <= boxsize; i++) {
	    temp = beta[1] + beta[2] * i - y[i];
	    chisq += temp * temp;
	}
    else
	for (i = 1; i <= boxsize; i++) {
	    temp = beta[1] + beta[2] * i - y[i];
	    for (j = 3; j <= nfit; j++) temp += beta[j] * x[i][j];
	    chisq += temp * temp;
	}
    return (chisq);
}

/* The functions below are based on those of the same names in Numerical
   Recipes (see above). */
void error(char error_text[])
{
    fprintf(stderr, "%s: %s\n", pname, error_text);
    exit(1);
}

double *vector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
    double *v = (double *)malloc((size_t)((nh-nl+2) * sizeof(double)));
    if (v == NULL) error("allocation failure in vector()");
    return (v-nl+1);
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
    int *v = (int *)malloc((size_t)((nh-nl+2) * sizeof(int)));
    if (v == NULL) error("allocation failure in ivector()");
    return (v-nl+1);
}

long *lvector(long nl, long nh)
/* allocate a long int vector with subscript range v[nl..nh] */
{
    long *v = (long *)malloc((size_t)((nh-nl+2) * sizeof(long)));
    if (v == NULL) error("allocation failure in lvector()");
    return (v-nl+1);
}

double **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
    long i, nrow = nrh-nrl+1, ncol = nch-ncl+1;
    double **m;

    /* allocate pointers to rows */
    m = (double **) malloc((size_t)((nrow+1) * sizeof(double*)));
    if (!m) error("allocation failure 1 in matrix()");
    m += 1;
    m -= nrl;

    /* allocate rows and set pointers to them */
    m[nrl] = (double *) malloc((size_t)((nrow*ncol+1) * sizeof(double)));
    if (!m[nrl]) error("allocation failure 2 in matrix()");
    m[nrl] += 1;
    m[nrl] -= ncl;

    for (i = nrl+1; i <= nrh; i++) m[i] = m[i-1]+ncol;

    /* return pointer to array of pointers to rows */
    return (m);
}

void free_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
    free(v+nl-1);
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
    free(v+nl-1);
}

void free_lvector(long *v, long nl, long nh)
/* free a long int vector allocated with lvector() */
{
    free(v+nl-1);
}

void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by matrix() */
{
    free(m[nrl]+ncl-1);
    free(m+nrl-1);
}