1 /* Copyright (C) 1997-1999 Adrian Trapletti
2
3 This library is free software; you can redistribute it and/or
4 modify it under the terms of the GNU Library General Public
5 License as published by the Free Software Foundation; either
6 version 2 of the License, or (at your option) any later version.
7
8 This library is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 Library General Public License for more details.
12
13 You should have received a copy of the GNU Library General Public
14 License along with this library; if not, write to the Free
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 time series bootstrap functions */
18
19
20 #include <math.h>
21 #include <R.h>
22
23
geodev(double p)24 static int geodev (double p)
25 /* Return geometric distributed random deviate with p(i) = p*(1-p)^i,
26 where 0<p<1 and i=0,1,..
27 See Fishman, G. S., 1996, Monte Carlo (Springer, NY, Berlin, Heidelberg), 221. */
28 {
29 double b, y, z;
30
31 b = (-1.0)/log(1.0-p);
32 y = exp_rand();
33 z = b*y;
34 return (int)z;
35 }
36
disuni(int n)37 static int disuni (int n)
38 /* Return discrete uniform distributed random deviate on {1,..,n} */
39 {
40 double temp;
41
42 temp = unif_rand()*(double)n+1.0;
43 return (int)temp;
44 }
45
WRAP(int i,int n)46 static int WRAP (int i, int n)
47 /* WRAP (i=..,-n+1,..,-1,0,1,..,n,n+1,.., n) = ..,1,..,n-1,n,1,..,n,1,.. */
48 {
49 if (i < 1)
50 return i%n+n;
51 else if (i > n)
52 return (i-1)%n+1;
53 else
54 return i;
55 }
56
StatBoot(double x[],double xBoot[],int n,double p)57 static void StatBoot (double x[], double xBoot[], int n, double p)
58 /* Generates a bootstrap sample xBoot[1..n] from x[1..n]
59 according to the stationary bootstrap resampling scheme
60 (Politis, D. N. and Romano, J. P., 1994, The Stationary Bootstrap,
61 J. Amer. Statist. Assoc. 89, 1303-1313).
62 Input is n, x[1..n], xBoot[1..n] and p, the parameter for
63 the geometric distribution, output xBoot[1..n]. */
64 {
65 int i, j, I, L;
66
67 i = 1;
68 while (i <= n)
69 {
70 I = disuni(n);
71 L = geodev(p);
72 j = 0;
73 while ((j < L) && (i <= n))
74 {
75 xBoot[i] = x[WRAP(I+j,n)];
76 i++;
77 j++;
78 }
79 }
80 }
81
BlockBoot(double x[],double xBoot[],int n,int L)82 static void BlockBoot (double x[], double xBoot[], int n, int L)
83 /* Generates a bootstrap sample xBoot[1..n] from x[1..n]
84 according to the blockwise bootstrap resampling scheme
85 (Kuensch, H. R., 1989, The Jackknife and the Bootstrap
86 for General Stationary Observations, Ann. Stat. 17, 1217-1241).
87 Input is n, x[1..n], xBoot[1..n] and L, the blocklength,
88 output xBoot[1..n]. */
89 {
90 int i, j, I;
91
92 i = 1;
93 while (i <= n)
94 {
95 I = disuni(n-L+1);
96 j = 0;
97 while ((j < L) && (i <= n))
98 {
99 xBoot[i] = x[I+j];
100 i++;
101 j++;
102 }
103 }
104 }
105
tseries_boot(double * x,double * xb,int * n,double * b,int * type)106 void tseries_boot (double *x, double *xb, int *n, double *b, int *type)
107 {
108 GetRNGstate();
109 if (*type == 0) StatBoot (x-1, xb-1, *n, *b);
110 else if (*type == 1) BlockBoot (x-1, xb-1, *n, (int)*b);
111 else error ("this type of bootstrap is not yet implemented\n");
112 PutRNGstate();
113 }
114
115
116
117