Estimates structural breaks in non-stationary univariate time series.
Y Array containing the time series. (Input)
NPCS Number of requested/estimated pieces or segments of the time series. NPCS is considered input when IFITONLY = 1. (Input/Output)
ARP A pointer to an array of size NPCS Χ 2. ARP is considered input when IFITONLY = 1. (Input/Output)
Column Index |
Description |
1 |
Requested/estimated break points |
2 |
AR order for each segment |
ARFIT A pointer
to an array of size NPCS Χ MAXARORDER containing
the AR coefficient estimates for each segment. ARFIT (i,j) is the j-th
coefficient for segment i where i=1, NPCS and j=l,
ARP(i,2).
(Output)
Note that the intercept is not reported.
ARSTAT A pointer to an array of size NPCS Χ 2. (Output)
Column Index |
Description |
1 |
Likelihood values for each of the fitted AR models |
2 |
Residual variances for each of the fitted AR models |
|
|
SCVAL Final value of the selection criterion. (Output)
MAXARORDERMaximum order to consider
for each AR model. (Input)
Default: MAXARORDER = 20.
IMTH Method of estimation. (Input)
Value |
Method | ||
0 |
Yule-Walker |
| |
1 |
Least Squares |
| |
2 |
Burg |
| |
Default: IMTH = 0.
ISELCRI Selection criterion. (Input)
Value |
Method |
0 |
Minimum Description Length (MDL) |
1 |
Akaikes Information Criterion (AIC) |
Default: ISELCRI = 0.
ILIKE Likelihood computation method. (Input)
Value |
Method |
0 |
Exact |
1 |
Approximate |
Default: ILIKE = 0.
IFITONLY Option to only fit the specified model. (Input)
Value |
Action |
0 |
No, do all the computations |
1 |
Yes, only fit the specified model |
Default: IFITONLY = 0.
IPRINT Printing option. (Input)
Value |
Action |
0 |
No printing |
1 |
Prints final results only |
2 |
Prints intermediate and final results |
Default: IPRINT = 0.
ISEED Seed of
the random number generator. (Input)
For the same data and parameter
settings, AUTO_PARM will return
the same results each time if ISEED = 1, 2
2147483646. If ISEED = 0, the system
clock will be used to generate a seed. The result will be
nondeterministic.
Default: ISEED =
0.
Note: The following arguments are for setting up and running the embedded Genetic Algorithm. In most situations, the default values should be used for these arguments. Users may wish to change some or all for testing or research purposes.
PDISTN Array of
length MAXARORDER + 1 giving
the probability distribution over the AR order variable p = 0,
,
MAXARORDER.
(Input)
j = 1,
, MAXARORDER + 1 is used
to randomly assign an AR order to breakpoint position j for a given
chromosome. PDISTN(j) >
= 0 and if SUM(PDISTN) is not equal
to 1, the values will be normalized, i.e., PDISTN(j)
= PDISTN(j)/SUM(PDISTN).
Default:
PDISTN
(j) = 1/(MAXARORDER + 1) for
all j.
MSPAN Array of
length MAXARORDER + 1
containing minimum number of observations required for valid estimates of AR
model with order p = 0,
, MAXARORDER.
(Input)
Default: MSPAN (p+ 1) =
2 *(number
of parameters) + 2 = 2 * (p + 2) + 2.
GAPARM Array of length 4 containing parameters that control the behavior of the genetic algorithm. These values should be strictly greater than zero and less than one to avoid unexpected results. (Input)
GAPARM (1)
Probability used to set initial break points in a chromosome.
Default:
MIN (MSPAN) / size(Y).
GAPARM (2)
Probability of crossover used to decide between a crossover and a
mutation.
Default: 1 MIN (MSPAN) / size(Y).
GAPARM (3)
In the mutation operation, probability an AR(p) model is enforced at the
current position.
Default: 0.4.
GAPARM (4)
In the mutation operation, probability a break point is disallowed at
the current position.
Default: 0.3.
Note: GAPARM(3) and GAPARM(4) must be valid probabilities and their sum must be between 0 and 1. 1 GAPARM(3) GAPARM(4) is the probability that the chromosome j inherits the parents chromosome j.
ISLAND Array of length 5 containing the migration policy parameters. (Input)
ISLAND (1)
Number of islands.
Default: 40.
ISLAND (2)
Number of generations that pass before migration occurs. Note
that the convergence of the algorithm is revised whether migrations take place
or not (see argument ISLAND(5)).
Default:
5.
ISLAND (3)
Number of subjects that migrate at each migration
event.
Default: 2.
ISLAND (4)
Population size (number of chromosomes) per island.
Default:
40.
ISLAND (5) =
Migration flag. If 1, migration is performed.
Default:
1.
MAXMIG Maximum
number of times that migrations may take place before the routine is stopped if
convergence has not occurred. (Input)
Default: MAXMIG = 20.
STOPITERS
Number of iterations. The routine will declare convergence and stop the
iterations if the criterion value (MDL/AIC) has not changed after STOPITERS consecutive
migrations. Otherwise, the algorithm will declare non-convergence and stop
after MAXMIG
migrations have taken place. See also MAXMIG and ISLAND(2). Note
that logically, STOPITERS < MAXMIG.
(Input)
Default: STOPITERS = 10.
Generic: CALL AUTO_PARM (Y, NPCS, ARP, ARFIT, ARSTAT, SCVAL [, ])
Specific: The specific interface names are S_AUTO_PARM and D_AUTO_PARM.
Routine AUTO_PARM estimates the structural breaks of a non-stationary time series using, with permission from the Authors, the method developed by Davis et al (2006). AUTO_PARM estimates a partition of the time index and models the time series in each segment as a separate auto-regressive (AR(p)) process. The routine returns the estimated breakpoints, the estimated AR(p) models, and supporting statistics.
For the observed time series ,
the problem is to find
m(m+1 = NPCS),
the number of breaks, their locations,
(ARP(:,1)),
and
(ARP(:,2)),
, the order of the AR
process in which the j-th segment is modeled. That is,
for
(for convenience,
and
) where
is an AR process of order
and , the noise sequence, is i.i.d.
with mean 0 and variance 1. Note that a series with m breaks will have
m + 1 segments.
The vector completely specifies a
piecewise AR model. To estimate this vector AUTO_PARM
optimizes, with respect to this vector, one of two selection critieria:
the first is a Minimum Description Length (MDL) criterion, and the second
is the Akaikes Information Criterion (AIC). The MDL is defined as
while the AIC criterion is given by
where, given a candidate value of the vector,
is the likelihood of
the fitted piecewise AR model evaluated at the parameter estimates,
The parameters of the j-th AR segment are estimated
by the choice of one of three estimation methods: Yule-Walker, Burg, or
Least Squares.
For simplicity, assume the mean of each series is 0 and that the errors are Gaussian. Then, the
piecewise AR model has Gaussian likelihood
where is the
variance-covariance of the j-th AR segment (of order
) and
is the vector of
observations of the j-th segment, i.e.,
To find the minimizer of either MDL or AIC, AUTO_PARM employs a
Genetic Algorithm with islands, migration, cross-over and mutations. See
Davis et.al. (2006) for further details.
AUTO_PARM approximates locally stationary time series by independent auto-regressive processes. Experimental results suggest that AUTO_PARM gives reasonable estimates of the structural breaks of a given time series, even if the segment series are not autoregressive. Also, based on experimental results, MDL gives better results than AIC as a selection criterion.
Informational error
Type Code
3 1 ISEED has been set out of range. ISEED is being reset to 123457.
3 2 MAXMIG migrations were reached in the genetic algorithm before the selection criterion value converged. Try increasing MAXMIG or using the double precision routine.
The examples below illustrate different scenarios using AUTO_PARM. The example series used in each case is the airline demand data (Box, Jenkins and Reinsel, 1994), which gives monthly total demand for the period January 1949 through December 1960. Each scenario sets the optional argument, ISEED = 123457.
use gdata_int
use auto_parm_int
use diff_int
use umach_int
implicit none
! Specifications for local variables
! arguments for auto_parm
integer :: npcs, maxarorder, ifitonly, iprint, island(5), iseed
integer, pointer :: arp(:,:)
real(kind(1e0)) :: x(144,1), scval
real(kind(1e0)), pointer, dimension(:,:) :: arfit, arstat
! Other locals
integer :: n, nvar, idata, nlost, iper(1), iord(1), nout
real(kind(1e0)) :: dx(144,1)
call umach(2, nout)
! Get the data
iseed = 123457
idata = 4
nullify(arp)
nullify(arfit)
nullify(arstat)
call gdata(idata, x, n, nvar)
! Example 1: Use defaults and print
! final results
write(nout,*)'Example 1: Use defaults'
write(nout,*)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
iseed=iseed, iprint=1)
! Example 2: Differenced series
! set period for the difference.
! iper is in years for this data set
write(nout,*)
write(nout,*)'Example 2: Differenced series'
write(nout,*)
! Set the order for the difference.
iper(1) = 1
iord(1) = 1
! Get differenced series dx
call diff(x(:,1),iper,iord,n,dx(:,1),nlost=nlost)
! Compare results on the differenced
! series
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(dx((1+nlost):n,1), npcs, arp, arfit, arstat, &
scval, iseed=iseed, iprint=1)
! Example 3: Original series
! lower maximum AR order
write(nout,*)
write(nout,*)'Example 3: Original series, lower order allowed'
write(nout,*)
maxarorder = 5
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
maxarorder=maxarorder, iseed=iseed, iprint=1)
! Example 4: differenced series, lower
! maximum AR order
write(nout,*)
write(nout,*)'Example 4: Differenced series, lower maximum'// &
' AR order'
write(nout,*)
deallocate(arp)
deallocate(arfit)
deallocate(arstat)
call auto_parm(dx((1+nlost):n,1), npcs, arp, arfit, arstat, &
scval, maxarorder=maxarorder, iseed=iseed, &
iprint=1)
! Example 5: Original series, force
! fit the segments
! Fit the specified model only at the
! break points
write(nout,*)
write(nout,*)'Example 5: Original series'
write(nout,*)'Force fit the segments at the break points'
write(nout,*)
npcs = 2
deallocate(arp)
allocate(arp(npcs,2))
arp(1,2) = 2
arp(2,2) = 1
arp(1,1) = 1
arp(2,1) = 60
deallocate(arfit)
deallocate(arstat)
call auto_parm(x(:,1), npcs, arp, arfit, arstat, scval, &
ifitonly=1, iseed=iseed, iprint=1)
end
Example 1: Use defaults
============== final results ===============
number of pieces: 2
selection criteria value: 684.242
total time: 0.7198125 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 1 0.77542
186.945 355.025
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
44 13 1.03701
-0.07802
-0.03890
-0.03452
0.11961
-0.12852
0.01990
-0.04886
0.08090
-0.13118
0.22122
0.53863
-0.61515
486.665 691.471
Example 2: Differenced series
============== final results ===============
number of pieces: 1
selection criteria value: 624.284
total time: 0.7204375 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 12 -0.02842
-0.22436
-0.16846
-0.24267
-0.10573
-0.22429
-0.12126
-0.26446
-0.07087
-0.24327
-0.07136
0.57129
619.321 297.351
Example 3: Original series, lower order allowed
============== final results ===============
number of pieces: 2
selection criteria value: 705.297
total time: 0.3136875 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 1 0.89533
270.393 333.564
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
63 2 1.19788
-0.35922
424.270 1632.337
Example 4: Differenced series, lower maximum AR order
============== final results ===============
number of pieces: 2
selection criteria value: 698.359
total time: 0.2981875 conv: 1
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 0 -------
335.565 357.388
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
77 1 0.33310
352.175 1786.345
Example 5: Original series
Force fit the segments at the break points
============== final results ===============
number of pieces: 2
selection criteria value: 712.521
==================== final model estimates =====================
break point order est. coeff. likelihood resid. var
ARP( 1,1) ARP( 1,2) ARFIT( 1,:) ARSTAT( 1,1) ARSTAT( 1,2)
1 2 1.12156
-0.24876
258.192 313.889
ARP( 2,1) ARP( 2,2) ARFIT( 2,:) ARSTAT( 2,1) ARSTAT( 2,2)
60 1 0.88605
443.696 1937.635
PHONE: 713.784.3131 FAX:713.781.9260 |