Bird's Zones
From BemlarGroup
[edit] Group
Peter Bird
Annie Chu
Dave Jackson
Yan Kagan
Rick Schoenberg
Qi Wang
Max Werner
Jiancang Zhuang
[edit] MODELS
| Model | Fast/Progr. | Boundary effect | Geometry | Focal Mech. | Magn. Scaling | Mag/time mask | Likl. Search |
|---|---|---|---|---|---|---|---|
| Agnes | Fast (F) | Boundaries | 2D-planar | No | S-Y, T-N | Triang. | Simplex |
| Max/Agnes | Fast (F) | Boundaries | 2D-planar | No | S-Y, T-N | Triang. | Simplex |
| Kagan/Jackson | Fast (F) | No Boundaries | 3D-spher. | Yes | S-Y, T-Y | Rectang. | N-R |
| Zhuang/Ogata | Fast | Boundaries | 2D-planar | No | S-Y, T-N | No | -- |
| Veen/Rick | Slow (R) | Boundaries | 2D-planar | No | S-Y, T-N | No | N-R |
| Annie/Rick | Slow (J) | Boundaries | 2D-spher. | No | S-Y, T-N | No | N-R |
| Console | Fast (F) | Boundaries | 2D-planar | No | S-Y, T-N | No | -- |
Fast/Progr.: F -- FORTRAN, R -- Statistical software, J -- Java; Magnitude Scaling: S-Y -- Space/Yes, T-N -- Time/No; Mag/time mask: A magnitude/time window after a strong earthquake to delete aftershocks that are incompletely sampled due to coda waves and overlapping records (Kagan, BSSA, 2004 [1]; Kagan and Houston, GJI, 2005 [2]; Helmstetter et al., BSSA, 2006 [3]; Ogata and Katsura, GRL, 2006); Likl. Search: Likelihood optimization technique, N-R -- Newton-Raphson (needs derivatives of likelihood function), Simplex (does not need derivatives).
[edit] Peter Bird
[edit] Bird's proposed tectonic zones
Defining tectonic zones
Here you can download data file for Bird's zone categories [4]
Format of PB2002_tectonic_zones.grd
[edit] Sample subcatalogs of CMT
Here are 5 subcatalog files derived from the shallow (depth <= 70 km below MSL)
and complete (m > 5.600; M > 2.818E17 N m = 2.818E24 dyne cm; note "16.05 conversion")
part of the Global CMT catalog, 1982.01.01-2008.03.31 inclusive.
All subcatalogs are in Bird's .eqc format, documented in an
Appendix
to Bird & Kagan [2004].
Zone 0 ( 256 EQs)
MAP
Zone 1 ( 870 EQs)
MAP
Zone 2 ( 500 EQs)
MAP
Zone 3 ( 706 EQs)
MAP
Zone 4 (4283 EQs)
MAP
[edit] Jiancang Zhuang
The differences between the ETAS models used by Zhuang and Ogata can be viewed within these 2 papers:
(1) Ogata Y. and Zhuang J. (2006). Space--time ETAS models and an improved extension. Tectonophysics. 413(1-2), 13-23. http://bemlar.ism.ac.jp/zhuang/pubs/ogata2006tectno.pdf
(2) Zhuang J., Chang C.-P., Ogata Y. and Chen Y.-I. (2005). A study on the background and clustering seismicity in the Taiwan region by using a point process model. Journal of Geophysical Research, 110, B05S18, doi:10.1029/2004JB003157. http://bemlar.ism.ac.jp/zhuang/pubs/zhuang2005jgr.pdf
(3) Ogata Y. (1988) Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes (in Applications). Journal of the American Statistical Association, 83 (401): 9-27. http://bemlar.ism.ac.jp/zhuang/Refs/Refs/ogata1988.pdf
(4) Ogata, Y., 1998. Space-time point-process models for earthquake occurrences, Ann. Inst. Statist. Math., 50(2), 379-402.
[edit] Yan Kagan
[edit] Aims of the project
I feel that we need to start by investigating how different branching models of seismicity approximate global earthquake catalogs. Global catalogs have an advantage of being considerably less inhomogeneous than local catalogs, and there are no spatial boundary effects which greatly complicate the analysis of local catalogs. Moreover, local seismicity is controlled by a few aftershock sequences of strong events, like the M7.5 1952 Kern County and the M7.3 1992 Landers earthquakes in southern California. On the other hand, it is important to analyze local catalogs as well, because, if we see that model parameter values are similar for worldwide and local catalogs (as was done in Kagan, 1991), this would present a good argument that the models are appropriate approximation of an earthquake occurrence and therefore they can be used for earthquake forecasting. Analogously, it is important to investigate various parametrizations of the branching models, especially the spatial and temporal fitting of seismicity patterns (see again Kagan, 1991) in order to find the best algorithms. Recently many researchers have applied branching models (usually various variants the ETAS model, see below) to approximate seismicity, but as far as I know, they have not worked with global catalogs. Apparently, only local catalogs in California, Japan, and Italy have been studied. Even those studies were conducted without any attempts to compare the results for different tectonic regions, various parametrization, and different catalogs.
[edit] A few remarks on the ETAS model of seismicity
In our joint project we would like to compare the results of fitting various branching models of earthquake occurrence to global catalogs. Therefore, I would like to mention a few problems that I see with the ETAS model. These problems may make such a comparison difficult.
1. The temporal influence for dependent shocks in the ETAS scheme is modeled by Utsu-Omori's dependence 1/(t + c). There are two problems caused by the use of such dependence. The first problem is that whereas the equation may be a reasonable approximation for fitting one aftershock sequence, many evaluations of the "c"-parameter show that it varies strongly, by orders of magnitude, even for mainshocks of similar size. Moreover, the c-estimates vary by an order of magnitude even for the same mainshock when using different data (catalogs). Kagan (2004) discusses the reasons for this phenomenon: the frequency range of a network, the distance to the closest station, installation of temporary stations following a strong event, workforce availability, etc. This effect may be of little importance in the consideration of individual sequences, but in the statistical analysis of earthquake catalogs we should expect that the parameters of the model have a property of statistical stability. Otherwise, derived parameter values are an average of various non-pertinent factors and as such may not be reproducible. In addition, many researchers interpret the c-value as having a physical significance; they should then be able to explain why "c" estimates are so unstable. This is in contrast to other seismicity parameters, like the p-value of Omori's law or the b-value of the Gutenberg-Richter relation, which hardly (if at all) change.
The other problem is that the coefficient "c" in the present ETAS implementations does not seem to be dependent on the magnitude of the preceding event. In reality, a relative temporal decrease of the aftershock numbers modeled by "c", varies strongly with the magnitude of a mainshock or a previous event. In my 1991 paper, I propose the following relation (see also Ogata, 1998, p.382 bottom):
tM = trM1 / 3 (8)
where tr is the coda duration time of an earthquake with the reference seismic moment Mr (tr = 3.45.10 − 3 days -- about 5 min -- for the magnitude 4 event). For a large earthquake this "dead" time would increase significantly: for M6 it is about 50 min, and for M8 it is eight hours.
If the c-value is independent of magnitude, then, for example, for a catalog with magnitude threshold 4.0, the c-estimate would be close to 5 min. This means that if the ETAS program is used for prediction, it would forecast after a M8 event a very large number of aftershocks in the few first minutes. It takes about one minute to rupture M8 fault, thus the forecast would be grossly off mark.
To avoid the bias, in the 1991 paper I remove from the catalog all the events within this time window, making appropriate changes in the likelihood function calculations. This procedure is similar to our removal of small earthquakes below the magnitude threshold from a catalog before its analysis.
It should be noted that Ogata (1983) introduced a coefficient "S" to describe a practical absence of aftershocks in the very initial part of a sequence. His estimate of S is 1/2 day (p. 120). Kagan (2004), Helmstetter et al. (2006), Ogata and Katsura (2006) suggest that the c-value does not only depend on the mainshock magnitude but on the size of the aftershock as well. For larger aftershocks the temporal delay is smaller. If these effects are not taken into account, the results of the statistical analysis of catalogs may be strongly biased.
2. As far as I can understand, the ETAS model presumes that a subsequent dependent earthquake may have a larger size than the previous event which caused it. In this case the name "Epidemic Type Aftershock Sequence" is misleading and it may be confusing, since the foreshocks are modeled as well.
REFERENCES:
Helmstetter, A., Y. Y. Kagan, and D. D. Jackson, 2006. Comparison of short-term and time-independent earthquake forecast models for southern California, Bull. Seismol. Soc. Amer., 96(1), 90-106. [5]
Kagan, Y. Y., 1991. Likelihood analysis of earthquake catalogues, Geophys. J. Int., 106, 135-148. [6]
Kagan, Y. Y., 2004. Short-term properties of earthquake catalogs and models of earthquake source, Bull. Seismol. Soc. Amer., 94(4), 1207-1228. [7]
Ogata, Y., 1983. Estimation of the parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure, Journal of Physics of the Earth, 31, 115-124.
Ogata, Y., 1998. Space-time point-process models for earthquake occurrences, Ann. Inst. Statist. Math., 50(2), 379-402. [8]
Ogata, Y., and K. Katsura, 2006. Immediate and updated forecasting of aftershock hazard, Geophys. Res. Lett., 33, L10305, doi:10.1029/2006GL025888.
[edit] References
Here are my references:
(1) Time-branching model (final research paper). Kagan, Y. Y., 1991. Likelihood analysis of earthquake catalogues, Geophys. J. Int., 106, 135-148. http://moho.ess.ucla.edu/~kagan/GJI_1991_like.pdf
(2) Time-branching model (first implementation). Kagan, Y. Y., and Knopoff, L., 1987. Statistical short-term earthquake prediction, Science, 236, 1563-1567. http://scec.ess.ucla.edu/~ykagan/sci1987_index.html
(3) Energy-branching model. Kagan, Y. Y., 1973. Statistical methods in the study of the seismic process (translated from Russian by D. Vere-Jones; with discussion: comments by M. S. Bartlett, A. G. Hawkes, and J. W. Tukey), Bull. Int. Statist. Inst., 45(3), 437-453 (scanned copy). http://moho.ess.ucla.edu/~kagan/Kagan_1973b.pdf
(4) Earthquake prediction using energy-branching model: Kagan, Y., and Knopoff, L., 1977. Earthquake risk prediction as a stochastic process, Phys. Earth Planet. Inter., 14, 97-108. http://moho.ess.ucla.edu/~kagan/pepi77.pdf
(5) See also discussion of these and other models in Kagan, Y. Y., 2006. Why does theoretical physics fail to explain and predict earthquake occurrence?, in: Modelling Critical and Catastrophic Phenomena in Geoscience: A Statistical Physics Approach, Lecture Notes in Physics, 705, pp. 303-359, P. Bhattacharyya and B. K. Chakrabarti (Eds), Springer Verlag, Berlin--Heidelberg; (in particular Fig. 17 and thereafter) http://scec.ess.ucla.edu/~ykagan/india_index.html
[edit] Creating subcatalogs
FORTRAN program to create smaller ARRAY with CHARACTER*1 variables: http://jumpy.igpp.ucla.edu/~kagan/PB2002_TECTONIC_ZONES1.FOR
FORTRAN program to create subcatalogs for each of Bird's tectonic categories: http://jumpy.igpp.ucla.edu/~kagan/BIRD_ZONES.FOR
[edit] Preprocessing catalogs
PDE catalog: FORTRAN program to correct for saturation of m_b and M_S magnitudes, take the average of four magnitudes (see Kagan, Y. Y., and Knopoff, L., 1980. Dependence of seismicity on depth, Bull. Seismol. Soc. Am., 70, 1811-1822) -- http://jumpy.igpp.ucla.edu/~kagan/STEP6.FOR
PDE and CMT catalogs: FORTRAN program to remove of close (in time) earthquakes from catalog (see Kagan, Y. Y., 1991. Likelihood analysis of earthquake catalogues, Geophys. J. Int., 106, 135-148; Kagan, Y. Y., 2004. Short-term properties of earthquake catalogs and models of earthquake source, Bull. Seismol. Soc. Amer., 94(4), 1207-1228) -- http://jumpy.igpp.ucla.edu/~kagan/STEP7A.FOR
P.S. In the preliminary results posted below for the PDE and CMT catalogs, these corrections have not yet been done. I hope that as the result of these corrections only minor changes in parameter values would occur.
[edit] Likelihood analysis of catalogs
FORTRAN program to calculate spatial moment of earthquake distribution in the PDE catalog, it is used later to normalize spatial kernel function (see Kagan, Y. Y., 1991. Likelihood analysis of earthquake catalogues, Geophys. J. Int., 106, 135-148): http://jumpy.igpp.ucla.edu/~kagan/PDEP2B.FOR
FORTRAN program to search for the maximum of the likelihood function, calculate values of parameters and information score for each of subcatalogs: http://jumpy.igpp.ucla.edu/~kagan/NLR1.FOR
Log file of FORTRAN program calculating spatial moment of earthquake distribution for the full PDE catalog M>=5: http://jumpy.igpp.ucla.edu/~kagan/PDEP2B.LOG
Log file of FORTRAN program searching for the maximum of the likelihood function, calculate values of parameters and information score for the full PDE catalog M>=5: http://jumpy.igpp.ucla.edu/~kagan/NLR1.LOG
Log file of FORTRAN program searching for the maximum of the likelihood function, calculate values of parameters and information score for the M>=5 PDE catalog in Subduction zones (Category 4): http://jumpy.igpp.ucla.edu/~kagan/NLR1_4.LOG
The FORTRAN programs for the CMT catalog and subcatalogs are similar.
A few technical remarks. Distances in global catalogs should be calculated as spherical lengths along great circles which requires more CPU time than simple computations for local catalogs. The events in global catalogs number in tens of thousands (see below), so calculations are extensive and the choice of appropriate software and optimization of programs should be carefully attended to. However, even on my old Alpha (200 4/233) computer that I've used for ten years, one likelihood iteration with 50000 events takes about 1/2 hour and about 100 iterations are needed to find a maximum of a likelihood function (see below). So this is doable. The software used is FORTRAN77 [9]. I spent some time 15-20 years ago optimizing the programs. If Moore's law is correct -- computer speed doubles every 2 years -- on modern computers one iteration should take about 1 min of CPU time.
[edit] CMT catalog
You can download CMT catalog in short format (I extracted information needed for likelihood function computation): http://jumpy.igpp.ucla.edu/~kagan/FPSHPR.DAT
Short version of the CMT catalog explained: http://jumpy.igpp.ucla.edu/~kagan/fpsh.txt
You can download CMT catalog in full (old and new) original formats: http://www.globalcmt.org/CMTfiles.html
To make it clear how I calculate the likelihood function, I have
posted several programs and input/output files I used to
process the CMT catalog. See my website
http://scec.ess.ucla.edu/~ykagan/like_index.html
The catalogs with calculated independence probability are stored at
http://scec.ess.ucla.edu/~ykagan/predictions_index.html
PARAMETER VALUES FOR VARIOUS SUBDIVISIONS OF CMT CATALOG,
1982--2007/03/31, Mw>=5.6
All 4 1 0 2 3
1. N 7720 5022 1004 279 584 831 2 Mmax 9.07 9.07 8.28 8.15 7.67 7.15 3. Inf/N 1.03 1.18 1.11 0.59 0.30 0.44 4. Ind/N 0.745 0.690 0.819 0.866 0.935 0.941 5. μ 0.131 0.169 0.093 0.099 0.042 0.035 6. b 1.02 0.98 0.96 1.05 1.28 1.23 7. c 0.64 0.59 0.56 0.45 0.52 0.59 8. θ 0.12 0.11 0.27 0.1* 0.27 *1.0 9. σ 0.30 0.29 0.15* 0.47 0.17 0.27 10. εr 21.8 22.1 21.6 18.0 19.5 17.7 11. εh 3.5 4.4 3.0* 5.1 3.0* 3.0*
This Table is similar in format to Table 3 in Kagan (Likelihood analysis of earthquake catalogues, 1991, GJI, 106, 135-148; http://scec.ess.ucla.edu/~ykagan/like_index.html). Briefly, N -- is the earthquake number, Mmax maximum magnitude in a subcatalog, Inf/N -- information score in bits/eq (1 bit means that the statistical uncertainty can be reduced on average by a factor of 2), Ind/N -- ratio of independent events to total, μ -- branching ratio (more branching - more dependent events), b -- b-value, c -- parent productivity exponent, δ = c / 1.5 (see Eq. 11), θ -- temporal exponent (1+θ is similar to "p" in Omori's law), σ -- focal size for M4 earthquake in km (sr in Kagan, 1991), εr -- horizontal error in km, εh -- vertical error in km (not informative for shallow earthquakes in the CMT catalog, since many events are assigned 10 or 15 km depth). The values preceded by an asterisk (*) mean that the parameter could not be evaluated and are assigned the upper limit, the values following by (*) are assigned the lower limit. Column numbers correspond to the tectonic categories (Bird's zones).
COMPARISON OF RESULTS FOR FULL CATALOG AND WITH CLOSE AFTERSHOCK REMOVED (3.3%)
All All/c 4 4/c
1. N 7720 7471 5022 4845 2 Mmax 9.07 9.07 9.07 9.07 3. Inf/N 1.03 0.857 1.18 0.997 4. Ind/N 0.745 0.758 0.690 0.707 5. \mu 0.131 0.131 0.169 0.168 6. b 1.02 1.01 0.98 0.98 7. c 0.64 0.62 0.59 0.57 8. \theta 0.12 0.1* 0.11 0.1* 9. \sigma 0.30 0.28 0.29 0.29 10. \eps_r 21.8 21.9 22.1 22.1 11. \eps_h 3.5 3.4 4.4 4.2
Close aftershocks have a strong influence on the information score, with about 3% aftershocks deleted, the score decreases by more than 15%.
[edit] PDE catalog, Maximum magnitude, M>=5.3
The catalog can be extracted from the NEIC web site http://neic.usgs.gov/neis/epic/epic_global.html . There are several formats available, some experimentation would be needed to find the appropriate one.
PARAMETER VALUES FOR VARIOUS SUBDIVISIONS OF PDE CATALOG,
1968--2007/01/01, M>=5.3
All 4 1 0 2 3
1. N 22515 14827 3583 960 1513 1632 2. Mmax 8.80 8.80 8.50 8.45 7.60 7.30 3. Inf/N 1.71 1.83 2.01 1.08 0.80 0.71 4. Ind/N 0.733 0.714 0.720 0.778 0.890 0.917 5. \mu 0.102 0.102 0.127 0.172 0.102 0.067 6. b 1.07 1.07 1.00 1.17 1.25 1.10 7. c 0.69 0.71 0.59 0.42 0.10 0.24 8. \theta 0.39 0.47 0.27 0.1* 0.54 *1.0 9. \sigma 0.37 0.36 0.30 0.26 0.15* 0.26 10. \eps_r 8.9 8.9 8.1 9.2 11.2 12.0 11. \eps_h 4.6 5.8 4.4 3.0* 3.0* 3.0*
This Table is similar in format to Table 3 of Kagan (Likelihood analysis of earthquake catalogues, 1991, GJI, 106, 135-148; http://scec.ess.ucla.edu/~ykagan/like_index.html). Briefly, M is the maximum magnitude of 4 given in the catalog for each earthquake (usually, m_b, M_S, M_L, and M_w), N -- is the earthquake number, Mmax maximum magnitude in a subcatalog, Inf/N -- information score in bits/eq (1 bit means that the statistical uncertainty can be reduced on average by a factor of 2), Ind/N -- ratio of independent events to total, \mu -- branching ratio (more branching - more dependent events), b -- b-value, c -- parent productivity exponent (see Eq. 11), \theta -- temporal exponent (1+\theta is similar to "p" in Omori's law), \sigma -- focal size for M4 earthquake, \eps_r -- horizontal error, \eps_h -- vertical error (not informative for shallow earthquakes in the PDE catalog, since many events are assigned 33 km depth). The values preceded by an asterisk (*) mean that the parameter could not be evaluated and are assigned the upper limit, the values following by (*) are assigned the lower limit. The column numbers correspond to the tectonic categories (Bird's zones).
[edit] PDE catalog, Maximum magnitude, M>=5.0
PARAMETER VALUES FOR VARIOUS SUBDIVISIONS OF PDE CATALOG,
1968--2007/01/01, M>=5.0
All 4 1 0 2 3
1. N 45942 29980 7686 2191 3296 2789 2. Mmax 8.80 8.80 8.50 8.45 7.60 7.30 3. Inf/N 1.90 1.99 2.16 1.40 1.28 0.98 4. Ind/N 0.680 0.661 0.687 0.716 0.816 0.869 5. \mu 0.141 0.140 0.133 0.234 0.182 0.109 6. b 1.05 1.04 1.05 1.17 1.17 0.93 7. c 0.63 0.65 0.64 0.35 0.02 0.18 8. \theta 0.28 0.33 0.23 0.12 0.40 0.55 9. \sigma 0.38 0.37 0.31 0.24 0.15* 0.20 10. \eps_r 9.5 9.7 7.9 9.8 10.1 12.6 11. \eps_h 3.0* 4.5 4.3 3.0* 3.0* 3.0*
This Table is similar in format to Table 3 in Kagan (Likelihood analysis of earthquake catalogues, 1991, GJI, 106, 135-148; http://scec.ess.ucla.edu/~ykagan/like_index.html). Briefly, M is the maximum magnitude of 4 given in the catalog for each earthquake (usually, m_b, M_S, M_L, and M_w), N -- is the earthquake number, Mmax maximum magnitude in a subcatalog, Inf/N -- information score in bits/eq (1 bit means that the statistical uncertainty can be reduced on average by a factor of 2), Ind/N -- ratio of independent events to total, \mu -- branching ratio (more branching - more dependent events), b -- b-value, c -- parent productivity exponent (see Eq. 11), \theta -- temporal exponent (1+\theta is similar to "p" in Omori's law), \sigma -- focal size for M4 earthquake, \eps_r -- horizontal error, \eps_h -- vertical error (not informative for shallow earthquakes in the PDE catalog, since many events are assigned 33 km depth). The values preceded by an asterisk (*) mean that the parameter could not be evaluated and are assigned the upper limit, the values following by (*) are assigned the lower limit. Column numbers correspond to the tectonic categories (Bird's zones).
COMPARISON OF RESULTS FOR FULL CATALOG AND WITH CLOSE AFTERSHOCK REMOVED (6.5%)
PARAMETER VALUES FOR VARIOUS SUBDIVISIONS OF PDECORR CATALOG,
1968--2007/01/01, M>=5.0
All 4 1 0 2 3
1. N 42925 27648 7205 2127 3217 2728 2. Mmax 8.80 8.80 8.50 8.45 7.60 7.30 3. Inf/N 1.25 1.24 1.50 1.17 1.04 0.76 4. Ind/N 0.695 0.674 0.709 0.736 0.831 0.881 5. \mu 0.171 0.177 0.146 0.241 0.185 0.121 6. b 1.04 1.03 1.04 1.17 1.17 0.93 7. c 0.57 0.56 0.61 0.27 0.0* 0.0* 8. \theta 0.11 0.13 0.10 0.1* 0.32 0.39 9. \sigma 0.37 0.35 0.32 0.15* 0.15* 0.15* 10. \eps_r 9.7 9.8 7.8 9.7 10.0 12.5 11. \eps_h 3.0* 4.8 4.7 3.0* 3.0* 3.0*
Similarly to the CMT catalog close aftershocks have a strong influence on the information score, with about 6% aftershocks deleted, the score decreases by more than 35%. Other parameters change is not as large: the number of independent events increases by about 2%, the branching coefficient (\mu) increases by about 20%.
[edit] PDE catalog, Maximum magnitude, M>=4.7
PARAMETER VALUES FOR PDE CATALOG,
1968--2007/01/01, M>=4.7
All 4 1 0 2 3
1. N 88920 56656 16235 4492 6730 4807 2. Mmax 8.80 8.80 8.50 8.45 7.60 7.30 3. Inf/N 2.01 2.052 2.159 1.526 1.786 1.374 4. Ind/N 0.637 0.614 0.660 0.679 0.753 0.813 5. \mu 0.177 0.184 0.152 0.240 0.225 0.158 6. b 1.00 0.97 1.06 1.09 1.10 0.87 7. c 0.58 0.57 0.64 0.41 0.12 0.15 8. \theta 0.22 0.23 0.18 0.11 0.40 0.45 9. \sigma 0.36 0.40 0.34 0.21 0.15* 0.15* 10. \eps_r 14.0 14.3 11.9 13.5 14.0 16.2 11. \eps_h 3.0* 3.0* 3.0* 3.0* 3.0* 3.0*
[edit] PDE catalog, m_b magnitude, m_b>=5.3
PARAMETER VALUES FOR VARIOUS SUBDIVISIONS OF PDE CATALOG,
1965--2007/01/01, m_b>=5.3
All 4 1 0 2 3
1. N 19318 13241 3216 869 1043 949 2. m_b (max) 7.30 7.30 7.30 6.90 6.50 6.40 3. Inf/N 1.66 1.75 1.93 1.00 0.71 0.59 4. Ind/N 0.711 0.689 0.712 0.802 0.909 0.929 5. \mu 0.143 0.159 0.124 0.141 0.074 0.057 6. b 1.44 1.43 1.36 1.38 1.69 1.87 7. c 0.90 0.86 0.99 0.63 0.37 0.45 8. \theta 0.25 0.26 0.21 0.1* 0.62 0.79 9. \sigma 0.77 0.80 0.54 0.46 0.15* 0.56 10. \eps_r 10.1 10.3 9.7 8.1 12.3 10.9 11. \eps_h 5.3 5.7 6.7 3.0* 3.0* 3.0*
This Table is similar in format to Table 3 in Kagan (Likelihood analysis of earthquake catalogues, 1991, GJI, 106, 135-148; http://scec.ess.ucla.edu/~ykagan/like_index.html). Briefly, N -- is the earthquake number, m_b (max) maximum m_b magnitude in a subcatalog, Inf/N -- information score in bits/eq (1 bit means that the statistical uncertainty can be reduced on average by a factor of 2), Ind/N -- ratio of independent events to total, \mu -- branching ratio (more branching - more dependent events), b -- b-value, c -- parent productivity exponent (see Eq. 11), \theta -- temporal exponent (1+\theta is similar to "p" in Omori's law), \sigma -- focal size for M4 earthquake, \eps_r -- horizontal error, \eps_h -- vertical error (not informative for shallow earthquakes in the PDE catalog, since many events are assigned 33 km depth). The values preceded by an asterisk (*) mean that the parameter could not be evaluated and are assigned the upper limit, the values following by (*) are assigned the lower limit. Column numbers correspond to the tectonic categories (Bird's zones).
[edit] Likelihood Analysis Results
The Tables below demonstrate that the estimates of alpha and mu are negatively correlated, other parameters look less correlated during the search. We need to obtain the Hessian matrix of the parameter estimates or at least plot two-dimensional plots of parameter estimates similar as was done in Bird and Kagan [2004]. Advantage of such plots is that if the relation is non-linear or parameter value needs to be restricted (because, for example, it goes into negative domain or to infinity, etc.) such plots give us a better information than the second-derivative matrix.
Generally, I feel that whatever is the value of parameters at the likelihood function maximum, they would provide a good prediction the effectiveness of which is approximately equal to the information score (Inf/N in Tables above). From the values of the log likelihood function difference one sees that the Akaike Info Criterion (AIC) should not work well for earthquake data -- we would need hundreds of additional parameters to make a dent in Inf/N value. The reason for it is likely that the AIC was developed for regularly varying variables, whereas seismicity patterns are governed by heavy-tailed, power-law distributions.
Speculations about parameters may, on the other hand, give us little information. The model is only an empirical fit to data (of questionable quality), so drawing any 'physical' conclusions is at best tentative and can be waste of time.
Table of likelihood maximum search for M>=5.0 PDE catalog (all earthquakes) is shown here: http://jumpy.igpp.ucla.edu/~kagan/all.txt The b-value does not vary during the search. See log file of FORTRAN program: http://jumpy.igpp.ucla.edu/~kagan/NLR1.LOG. Postscript plots of F, mu, and alpha are shown in http://jumpy.igpp.ucla.edu/~kagan/PDE_5_ALL1.PS and http://jumpy.igpp.ucla.edu/~kagan/PDE_5_ALL2.PS
Table of likelihood maximum search for M>=5.0 PDE catalog in subduction zones is shown here: http://jumpy.igpp.ucla.edu/~kagan/subd.txt. b-value does not vary during the search. See log file of FORTRAN program: http://jumpy.igpp.ucla.edu/~kagan/NLR1_4.LOG. Postscript plots of F, mu, and alpha are shown in http://jumpy.igpp.ucla.edu/~kagan/PDE_5_SUB1.PS and http://jumpy.igpp.ucla.edu/~kagan/PDE_5_SUB2.PS
See our ms discussing the likelihood search results: http://scec.ess.ucla.edu/~ykagan/globe_index.html
[edit] Max Werner
This document contains a brief definition of the ETAS model, along with a derivation of the branching ratio from a pure and a truncated Gutenberg-Richter law. The branching ratio depends on the corner magnitude if calculated correctly. The file can be downloaded from [10].
Here is a copy of Ogata's 1998 paper:
Y. Ogata, Space-Time Point Process Models for Earthquake Occurrences, Ann. Inst. Statist. Math. Vol. 50, No. 2, 379-402, 1998.
[11]
[edit] Annie
[edit] Model
The ETAS model used here is from Ogata 1998, equation (2.4)
[edit] Data
Currently (as of 2007/11/09) I'm fitting the CMT data. I'm using dataset from Kagan: http://jumpy.igpp.ucla.edu/~kagan/FPSHPR.DAT
As for testing, I use a portion of the above dataset, for two reasons: 1. faster to run the program; 2. to get a idea on what initial values and range of possible values may work better without taking too much time running the large dataset. (Range of possible values are for the parameter a, p, and q, where a is alpha in the model's triggering function of the ETAS model; p is 1 + omega; q is 1 + rho.
[edit] Computation & Algorithm
The computational algorithm is an EM-type estimation procedure developed by Veen and Schoenberg (2007). Great circle distance is used within the spatial region: longitude = -180 to 180; latitude = -90 to 90.
Code is originally developed in R. I rewrote it in Java and that's the one I'm running.
[edit] Preliminary Results and Implementation Details
2008/05/17 Implementation of Spherical Polygon Area Computation
2008/05/04 Inhomogeneous Test on a region around China
2008/04/25 Homogeneous Test on several regions (Indonesia , Japan Honshu, Southern California, Taiwan)
2007/11/09 Homogeneous Test on several small sets
[edit] World Map + CMT + Bird's Zones
World Map plotted with zones only
World Map plotted with zones and CMT events (1976 to 2007)
[edit] Model Fitting of Global Data(CMT)
Results of fitting CMT data on Ogata's 1998 space-time model (2.4) conditional intensity lambda = background rate mu + summation of triggering function g where g = k * (t + c)^-p * ( R*R / exp(a * (M_I – M0) + d)^-q
Using time window = 1982/01/01 to 2007/03/31 Minimum magnitude cutoff = 5.55 Link to results: http://www.stat.ucla.edu/~chea/etas/CMTResult_1982TO2007_555_70.htm
[edit] Model Fitting of Global Data(PDE)
Results of fitting PDE data on Ogata's 1998 space-time model (2.4), using the same model as described above.
Time window: 1973/01/01 to 2006/12/31 Only shallow events (0 KM to 70 KM in depth) were used. Lower magnitude cutoff M0 = 5.0 Reference Model: stationary, homogeneous, Poisson
Link to results: http://www.stat.ucla.edu/~chea/etas/PDEResult_1973TO2006_50_70B.htm and plots based on the parameter estimates: http://www.stat.ucla.edu/~chea/etas/PDEResult_1973TO2006_Plot.htm
[edit] Poster presentation for StatSei6
Link to poster (6.9MB PDE file) http://www.stat.ucla.edu/~chea/etas/StatSeiPosterG.pdf
[edit] Tables
Table of MLE only:
http://www.stat.ucla.edu/~chea/global/MLE_Table_1973TO2006_50.doc
Table of summary or area, number of events, and table of MLE with SE:
http://www.stat.ucla.edu/~chea/global/Tables.pdf
[edit] Plots and figures
Plots to show distance containing 95% of the aftershocks vs. magnitude:
http://www.stat.ucla.edu/~chea/global/DistanceMag_123405.pdf (all of the 6 zones)
http://www.stat.ucla.edu/~chea/global/DistanceMag_12340.pdf (zones 1, 2, 3, 4 and 0)
Plots of MLE vs. zones, and pairs of MLE are under:
http://www.stat.ucla.edu/~chea/global/varplot/
A collection of figures with maps and some diagnostics plots (including estimated conditional intensity (in the space dimension), estimated g on R and R^2, estimated g. vs. empirical rate of aftershocks, isotropic triggering of aftershocks, etc.)
http://www.stat.ucla.edu/~chea/global/Figures.pdf
[edit] Qi
[edit] Sample subcatalogs of CMT
Here are 5 subcatalog files derived from the Global CMT catalog with condition: Time:1982.01.01-2007.03.31;Depth <= 70 km;Moment magnitude >5.6 or Moment >2.5119e+024.
Zone 0 ( 237 EQs)
Zone 1 ( 898 EQs)
Zone 2 ( 487 EQs)
Zone 3 ( 723 EQs)
Zone 4 ( 4407 EQs)
MAP
[edit] ETAS results of CMT
[edit] Acknowledgements
This work is supported by the Geophysics Program of the National Science Foundation under grant EAR-0711515 to the University of California. All opinions expressed are those of individual members of the research group, and do not relect official positions of the National Science Foundation, nor of the U.S. Government.
