Regression modeling method of space weather prediction

0.77 MB·English

by  Aleksei Parnowski

AstrophysicsandSpaceScience
DOI10.1007/s•••••-•••-••••-•
Regression modeling method of space weather prediction
A.S. Parnowski
0
1
0
2

n
a
J

2
1 (cid:13)c Springer-Verlag••••

] Abstract A regression modeling method of space like the geomagnetic indices. Since this article centers
h
p weather prediction is proposed. It allows forecasting on the variations of the geomagnetic field, the latter
- Dstindex upto 6hoursaheadwithabout90%correla- meaning will be used. The second part of this defini-
e
c tion. It can also be used for constructing phenomeno- tionindicatespracticalmanifestationsofspaceweather.
a logicalmodelsofinteractionbetweenthesolarwindand The impact of the space weather on technological sys-
p
the magnetosphere. With its help two new geoeffective tems is generally accepted (see Marubashi (1989)) due
s
. parameters were found: latitudinal and longitudinal to a number of spectacular events like the superstorm
s
c flow angles of the solar wind. It was shown that Dst of 1989 when Canada’s power grid was disabled for 9
si index remembers its previous values for 2000 hours. hours and numerous spacecraft failures due to “killer
y electrons” causing arcing in electronic components, see
h Keywords space weather; prediction; forecasting;
Romanovaet al. (2005). The impact on human health,
p magnetic storms; statistics; regression
[ however, is disputed by most specialists. Nevertheless,

thelatestreports(e.g. Khabarova&Dimitrova(2008),
2
v 1 Introduction Stoupel et al. (2006)) indicate that there is indeed a
1 strong correlation between the rate of sudden cardiac
7 The humankind studies space weather for more than deaths and the space weather.
2 4000yearsstartingfromthefirstmentionsofaurorasin The space weather problem is twofold. The first
3
ancient Chinese literature. The term “space weather” aspect is purely practical and aims for prediction
.
6 itself exists for almost a century. The official defini- and, eventually, mitigation of adverse effects of space
0
tion adopted by COSPAR states that “Space weather weather. Ideally, this task should be accomplished
9
0 describes the physical processes induced by solar ac- by launching a vast number of spacecraft which will
: tivity that have impact on our terrestrial and space monitorthe Sun-Earthregionforlarge-scalestructures
v
i environment, on ground based and space technological like CMEs. Unfortunately, the resources of the hu-
X systems, and on human activities and health.” The mankind are insufficient to produce and maintain such
r first part of this definition actually covers two spatial a large space fleet as well as to process all the data
a
scales of space weather, because when we speak about delivered by these spacecraft. So, today we should
space weather in space, e.g. in connection with space-
use the resources at hand, which include a few solar
craft failures, we usually mean some local parameters
wind spacecraft(ACE,WIND, SOHO, and STEREO),
of the environment, and when we speak about space
magnetosphericspacecraft(CLUSTER,THEMIS),and
weather on the Earth, e.g. in connection with human
ground-based stations (Intermagnet, MAGDAS, etc.),
health, we usually mean some integral characteristics
to develop forecast techniques that will be used in fu-
ture. Thus, we should try to predict space weather
A.S.Parnowski with what data we have, and we should aim for pos-
SpaceResearchInstituteofNASUandNSAU sibly longer prediction times to allow for some kind of
prosp. Akad. Glushkova, 40, korp. 4/1, 03680 MSP, Kyiv-187, countermeasures.
Ukraine
The second aspect is mostly academic and involves
tel: [phone], fax: [phone]
e-mail:[email protected] study of the processes in the near-Earth space and,
specifically, understanding of interaction between the
2
solar wind and the magnetosphere. Naturally, improv- 0.3
ACF
ing our knowledge of the underlying physics signifi-
cantly improves predictive capabilities, so fulfilling the
0.2
second task will significantly help with the first one.
Modern conceptions of solar wind-magnetosphere in-
teraction are mostly based on phenomenological mod-
0.1
elsconstructedin1960’s. However,therearenumerous
problems these models cannot answer. This is largely
due to the fact that these models were developed at 0
the very beginning of the space era when data quality 0 4000 8000 12000 16000 20000
andquantitywereimmeasurablyworsethantoday. For Time offset, hrs
morethan40yearswecollectedastonishingamountsof
Fig. 1 Autocorrelation function of Dst. Horizontal lines
data aboutsolarwindparametersandgeomagneticac-
correspond to top and mean incidental correlation levels in
tivity and now it is time to put them to good use.
abscence of periodic variations. Gray sine has a period of
1/2 year and depicts seasonal variations
2 Possible approaches to space weather
This method provides the shortest prediction time
prediction
(up to 1 hour) with moderate accuracy (∼ 70%). Po-
tentially this method could demonstrate far better re-
Spaceweatherpredictionisachallengingandnontrivial
sultsifthephysicsbehindthemagneticstormswasless
activity, see Li et al. (2003). The most straightforward
of a mystery.
approach to space weather prediction is studying the
4. Numerical modeling, see e.g. Liu et al. (2008),
whole complex chain of physical processes involved in
McKenna-Lawloret al. (2008).
magnetospheric dynamics and conjugating them in a
This method providesa goodpredictiontime (up to
globalmodeloftheevolutionofthemagnetosphereun-
severaldays)butitsaccuracyvariesinhugelimits. The
der the influence of the solarwind. Unfortunately, this
accuracyofthesemethodsislimitedbytheirinabilityto
is not yet possible due to our poor understanding of
correctlydescribeplasmainstabilities. Besides,thering
the physics of the interaction between the solar wind
current can not be described in the framework of ideal
and the magnetosphere. For this reason, different ap-
MHD,whichformsthebasisofmostnumericalmodels.
proaches should be tried.
However, they can adequately describe the motion of
AccordingtoKhabarova(2007),todaytherearesev-
e.g. magneticcloudsintheinterplanetaryenvironment,
eral established methods of space weather prediction,
but rely on different methods to detect them.
listed below.
5. Multidimensional time series analysis.
1. Morphologicalanalysis of solar images.
This method provides a moderate prediction time
This method provides the longest prediction time
(up to several hours) with the highest accuracy (>
(up to a week). Its accuracy is unknown since it is
80%). They are very effective and easy to use but
usedfortheacademicpurposesonly. Todayitispurely
stronglydependonsatellitedataavailability. Theseare
manual and thus almost useless for practical implica-
“blackbox”or“input-output”models,whichseekonly
tions.
toreproducethesystem’soutputinresponsetochanges
2. Detectionoflarge-scaleperturbationsinthe solar
of its inputs. The model terms are usually physically
wind, see e.g. Eselevich & Fainshtein (1993), Eselevich
interpretable and thus useful for construction of new
et al. (2009).
phenomenologicalmodels. Forthisreason,thismethod
This method provides a very good prediction time
cannotonlyprovidespaceweatherforecastper se,but
(up to several days), but is capable of predicting less
alsocanimproveourknowledgeofthe physicsinvolved
than 10% of the most intense storms. While it is very
and thus increase the efficiency of other methods.
inaccuratewhenusedalone,itcanprovetobeusefulin
Further we shallspeak aboutthe lastmethod, keep-
combinationwithoneofthefollowingshort-termmeth-
ing in mind that its results can be used later to assist
ods.
other methods. First of all, let us discuss its existing
3. Constructionofempiricalmodels,seee.g. Burton
implementations.
etal.(1975),Valdiviaetal.(1996),O’Brien&McPher-
Multidimensional time series analysis can be per-
ron(2000a),O’Brien&McPherron(2000b),Temerin&
formed using the methods of statistics, signal process-
Li (2002), Temerin & Li (2006), Ballatore & Gonzalez
ing, informatics, fuzzy logic etc. The most widely used
(2003), Cid et al. (2005), Siscoe et al. (2005).
variationsare artificial neuralnetworks(e.g. Kugblenu
Regressionmodelingmethod 3
600 30
Monte-Carlo simulations F
Normal distribution
400 20
s
nt
oi
p
of
er
b
m
u 99.95%
N
99.9%
200 10
99.5%
99%
97.5%
95%
90%
0 0
0 0.004 0.008 0.012 0.016 0 8 6 4 2 0 8 6 4 2 0 8 6 4 2 0 8 6 4 2
4 9 4 9 4 8 3 8 3 8 2 7 2 7 2 6 1 6 1
Correlation coefficient 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9
Time offset, hrs
Fig.2 DistributionofcorrelationcoefficientofDstatvery
Fig. 3 Dependence of Fisher significance F of the corre-
large time offsets in abscence of periodic variations spondingterm in equation (1) on the time offset for the1h
autocorrelation model
etal.(1999),Watanabeetal.(2002),Wingetal.(2005),
Pallocchia et al. (2006)), optimization (e.g. Zhou & of 5 intense CME driven storms during the 1996-2002
Wei (1998), Balikhin et al. (2001), Harrison & Drezet period using another 46 CME driven storms to train
(2001)), and correlation analysis (e.g. Rangarajan & his model.
Barreto (1999), Oh & Yi (2004), Wei et al. (2004), Hereafter we propose a new approach, named “re-
Johnson & Wing (2004), Johnson & Wing (2005)). gressionmodeling”, which already allows achieving ac-
Neural network approach provides short-term predic- curate (∼ 90%) 6 hours ahead forecast of the Dst in-
tions up to 4 hours with the correlation coefficient of dex, which we will use as a quantitative characteristics
0.79 in the paper by Wing et al. (2005). Earlier im- of space weather. This method can be easily extended
plementations of this approach experienced significant to predict other geomagnetic indices like Kp or Ap.
difficulties predicting strong geomagnetic storms with
Kp > 5, but this approach remains one of the most
popular alongside the empirical methods. Optimiza- 3 Description of the regression modeling
tion approach seems to be more successful being able method
to provide 8-hour predictions in the paper by Harrison
& Drezet (2001). However, in the papers based upon The proposed method is statistical, but has some fea-
the optimization methods the volume of the dataset tures of empirical models. It is based upon the regres-
used is insufficient to correctly describe secular varia- sion analysis and the mathematical statistics. In its
tionsofgeomagneticactivity. Correlationanalysisgives frameworkthepredictedDstvalueissoughtintheform
interesting results, but it was used solely for develop-
ing and constraining empirical models (see Johnson & Dst(j+k)=XCixi(j), (1)
Wing (2004)). However, most of these methods have i
a common feature: they lead to a regression relation-
wherej isthenumberofcurrentstep(numberofhours
ship at some point, so it seems natural to skip all the
since Jan 1, 1963), k is the prediction length, C are
i
preliminarystepsandinstantlyusetheregressionanal-
the regression coefficients, and x are the regressors,
i
ysiswithoutunnecessarymultiplicationofentities. Re-
which are functions and combinations of input quan-
gressionanalysisitself wasattempted earlierby Srivas-
tities, which are already measured at the time when
tava (2005), but it was used to estimate the probabil-
prediction is made. Values of C are determined by
i
ityofintense/super-intensestormoccurencedepending
leastsquaremethod(LSM)overalargesampleofsolar
on the solar and interplanetary parameters. Srivastava
wind and geomagnetic data (see next chapter), with
(2005) was able to predict 2 of 4 super-intense and 5
equal statistical weights of all points. The statistical
4
0 0.6
-5 T
st, nT Dst, n0.4
Mean D-10 nput in 0.2
-15 I
-20 0
0 30 60 90 120 150 180 210 240 270 300 330 360 0 30 60 90 120 150 180 210 240 270 300 330 360
Days of year Days of year
Fig. 4 Seasonal variation of Dst Fig. 6 Sumoftermsdirectlydescribingseasonalvariation
of Dst
-10
0.4
n Dst, nT-11 Dst, nT
Mea-12 ut in 0.2
p
n
I
-13
0 4 8 12 16 20 24 0
UT hours
0 4 8 12 16 20 24
UT hours
Fig. 5 Diurnalvariation of Dst
Fig. 7 Sum of terms directly describing diurnal variation
of Dst
significance of the regressorswas determined by Fisher
test (F-test) (see Fisher (1954), Hudson (1964)). This
ysis: it allows checking the statistical significance of
test allows separating significant and insignificant re-
any regressor, which can be useful for verifying differ-
gressors. Theinsignificantparametersarethenrejected
ent physical hypotheses. In this sense we will call a
andtheroutineisrepeateduntiltheregressioncontains
parameter “geoeffective” if it appears at least in one
onlysignificantregressors. Ofcourse,thismethoddoes
statistically significant regressor.
notguaranteethatallthe significantregressorswillen-
More details of this method can be found in the ar-
tertheregression,butphysicalconsiderationsandbrute
ticle Parnowski(2009a).
force in the form of trial and error provide us with re-
questedreliability. Theregressorsx aregenerallynon-
i
linear, so from the control theory’s point of view, this
4 Data and routine
method is able to describe discrete dynamical systems
with strong nonlinearity. This is an essentialfeature of
The OMNI2 (2009) database was used. It contains
the regressionmodeling method.
IMF, solar wind and geomagnetic data, averaged over
There is only one manual operation in this method
1-hour intervals (49 parameters in total, starting from
– selectionof regressorsto be considered. For this pur-
Jan 1, 1963). It was supplemented with provisional
pose all known models, basic physical considerations,
Dst data, taken from WDC for Geomagnetism (Ky-
and random choice are used. Naturally, common sense
oto). Thus a continuous 44-year Dst time series was
also counts: for example it would be silly to add IMF
obtained.
components in GSE and GSM coordinates at the same
We estimatedthe geoeffectivenessofaparameterby
time. If some regressorsx appeared to be statistically
i coefficientsandstatisticalsignificancesofallregressors,
significant,wealsocheckedthe significanceofproducts
which contain this parameter. This was done in the
of their powers xpi, where p can be any real num-
Q i i following way. After processing the data with the least
i
ber, including zero, but for practical purposes we used square method, Fisher significance parameter F was
integer values of pi in the range from 0 to 6. This determined for each regressor. All the F values were
yields a very important feature of the regression anal- comparedto the values 2.7055,3.84,5.02,6.635,7.879,
Regressionmodelingmethod 5
10.83 and 12.1, which correspond to statistical signif-
icance of 90, 95, 97.5, 99, 99.5, 99.9 and 99.95% re-
spectively. Then, insignificant regressors were rejected
and the routine was repeated until all the regressors
were significant. The number of significant regressors
depends on the selected significance threshold. All re-
sultsgivenhereincorrespondtothesignificancethresh-
old of 90%. In contrast to empirical models we do not
addfittingparametersandalltheregressorshavephys-
icalmeaning. Thedescribedroutinewasappliedtothe
sample,obtainedfromtheinitialdatasetafterrejecting
filled values. This sample can be divided into two sub-
samples, corresponding to quiet (Dst > −50nT) and
perturbed (Dst≤−50nT)conditions.
First, we determined which previous Dst values are
statistically significant. For this purpose, we con-
structed an autoregression (see details in Parnowski
(2009b))
N
Fig.8 TemporalvariationofDst. Darkerspotscorrespond
Dst(j+k)=C0+XCiDst(j−i+1), (2)
tolower values
i=1
where N is the “age” of the oldest Dst value. This
temporal variations. The ACF crosses the top line for
model alone is not sufficient to correctly predict Dst,
the first time at ∼ 6000 hours, though the difference
but it sets a basis for the construction of models that
between the ACF and the sine with a half-year period
are able to do so. Let us determine the maximum rea-
crosses it at ∼ 2000 hours, which is about 3.5 27-day
sonablevalueofN. Forthis purpose,we plotthe auto-
periods, so we will assume the latter value as a rough
correlation function (ACF) of the Dst index for k = 1
estimation of N. This hints that rather old Dst values
(see Fig. 1). One can see that ACF tends to a sinu-
canbequitesignificant. Besidesthehalf-a-yearperiod-
soid with a period close to half a year. This is caused
icity one can also notice the 27-day periodicity, caused
by seasonal variations. This yields a question: if there
by Carringtonrotation of the Sun, which can be taken
were no temporal variations, what would ACF tend to
into account by adding the sunspot number R to the
at large offsets? If the distribution of Dst was normal,
regression.
the answer would be zero. However,the distribution is
Let us return to equation (2). Applying the F-test
not normal, so ACF can tend to some non-zero quan-
we can determine which previous Dst values are sta-
tity. To determine this quantity we need to remove
tistically significant (see Fig. 3). We did not search
temporal variations. For this purpose we need to cal-
statistically significant Dst values for N > 900, but it
culate the ACFofa randomsample withthe samesta-
ispossiblethatthereareevenolderstatisticallysignifi-
tistical characteristics as the Dst sample. The easiest
cant values. A similar situation was reported by John-
way to get such a sample is to process the Dst sample
son & Wing (2004) regarding Kp: ”the significance is
with a permutation method, which is widely used for
often quite large for extended periods of time (10-20
determination of correlation functions, e.g. in astron-
days)”. As our analysis shows, Dst index can “remem-
omy. This method is based on random shuffle of the
ber” its previous values for significantly longer periods
sample. Usingthismethodmanytimes(10000timesin
of time. In fact, after adding regressors, correspond-
ourcase)andcalculatingthecorrelationcoefficienteach
ing to satellite data, some of the previous Dst values
time, we get the distribution of the correlation coeffi-
become insignificant. We found that after the addition
cient by Monte-Carlo method. The distribution of the
oftheseregressorstherearestillstatisticallysignificant
correlationcoefficientforthissample(Fig. 2)appeared
to be very close to a normal distribution with mean valuesasfaras801hoursago(33daysand9hours)for
0.008 and variance 5.1·10−6. The maximum recorded k = 1. The statistical significance of this oldest value
value in 10000 trials was equal to 0.015. The top and is over 99.9%.
the mean values are depicted on Fig. 1 by horizontal At this point we already have a large number of
lines. As one can see on Fig. 1, in reality the correla- regressors, describing just the previous Dst values
tion coefficient exceeds this value at most times due to (autoregression), without satellite data and nonlinear
6
30000 1
Input in Dst, nT
s
point20000 -15 -10 -5 5 10 (cid:84)V, °15
of
er -1
b
m10000
u
N
-2
00
-3
-5
T
st, n Fig.10 Sumoftermsdescribingthelatitudinalflowangle
D -10
n
a
e
M
-15
-20
-15 -10 -5 0 5 10 15
(cid:84), °
V
Fig. 9 Distribution of the latitudinal flow angle and the
corresponding mean Dst values
terms. If we add those, the number of regressors will
only increase.
AfterdeterminingwhichpreviousDstvaluesaresta- Fig. 11 Seasonal dependence of latitudinal flow angle’s
inputin Dst
tistically significant, we added all the solar wind pa-
rameters available in the OMNI 2 database. Then, we
added nonlinear terms as discussed in Section 3. Af- that this mechanism would give only 17% of observed
ter adding a new regressor,all the significances are re- asymmetry. Takalo & Mursula (2001) connected the
calculated, and some of the old regressors can become diurnal variations of Dst with an uneven distribution
insignificant. The total number of regressors is about of Dst network stations. Let us use this known effect
150-200. Since it is very large, we will not give here to validate our method.
any lists of regressors or coefficients even for the sim- If we select two subsamples, corresponding to sum-
plestcasek =1,butthepreliminarylistisgiveninthe mer and winter in northern hemisphere, bounded by
paper Parnowski(2008). vernalandautumnalequinoctia,andverifythehypoth-
esis that the difference between the corresponding av-
erage Dst values is statistically significant using a one-
5 Identification of new geoeffective parameters sidedStudenttest,weobtaint∞ =80.264,whichiswell
over 99.95% significant. Values of t∞ corresponding
In this section we will demonstrate how this method to 99 and 99.95% significance levels are equal to 2.334
canbeusedforidentificationofgeoeffectiveparameters. and 3.31 respectively. For diurnal asymmetry Student
We will use four parametersas anexample: DOY (day testgivest∞ =8.774,whichcorrespondstosignificance
of the year), UT (universal time) and latitudinal and levelof more than 99.95%. Note that formally Student
longitudinal flow angles of the solar wind. test is applicable when Dst is normally distributed. In
OnFig. 1 one couldsee a clearseasonaldependence fact, the distribution ofDst is slightly asymmetric,but
of the Dst index. This dependence was described in taking into account the obtained values of t∞, we can
many articles, for example by O’Brien & McPherron be sure in qualitative conclusions made. Figs. 4 and 5
(2002),Lyatskyetal.(2001),Takalo&Mursula(2001), show the histograms of seasonaland diurnal variations
andCliveretal.(2000),butthereasonbehinditisstill of Dst index.
disputed. Most authors believe these asymmetries are Taking this known geoeffective factor as an exam-
causedby either oftwo cusps turning to the sunlit side ple we demonstrate how easily one can take it into ac-
due to annual rotation of the Earth with respect to count using regression approach. To do so one should
the Sun. However, O’Brien & McPherron (2002) state simply add terms a1(j) = sin((j −1920)π/4383) and
Regressionmodelingmethod 7
40000 2
Input in Dst, nT
s
nt
poi 1
er of 20000
b
m
u
N
-15 -10 -5 5 10 15
(cid:77), °
00 V
-1
T
n Fig. 13 Sumoftermsdescribingthelongitudinalflowan-
Dst, -20 gle
n
ea Bin size = 1°
M
Quiet conditions (Dst > -50 nT)
All data
-40
-15 -10 -5 0 5 10 15
(cid:77), °
V
Fig.12 Distributionofthelongitudinalflowangleandthe
corresponding mean Dst values. Grey columns correspond
to quiet conditions, white columns – to all data
a2(j) = cos((j − 1920)π/4383) into the regression.
Here j is once again the number of hours since Jan
Fig. 14 Seasonal dependence of longitudinal flow angle’s
1, 1963, 1920 is the number of hours between the be-
inputin Dst
ginning of the year and the vernal equinox, and 4383
is the number of hours in half a year. The first of
these terms is significantanddescribes summer/winter term in brackets is of order (ω∆t)2, and the second –
asymmetry,andthesecondone(whichappearsstatisti- ω∆t in the natural assumption that ω∆t ≪ 1. So, it
cally insignificant) describes an absent spring/autumn will seem that the coefficient is Aω∆t rather than A.
asymmetry. Likewise, for diurnal asymmetry the cor- Notethatthisisjustanexampleandhasnothingtodo
responding terms will be b1(j) = sin((j −2)π/12) and with actual regressors.
b2(j) =cos((j−2)π/12). Here 2 is the time difference However if we look on the distribution of mean Dst
between UT and the northern geomagnetic pole, and values vs. time of the year (Fig. 4), we see a much
12 is the number of hours in half a day. Both these more complicated pattern of seasonal variations of ge-
terms are significant. omagnetic activity. Among other features there is a
Thecoefficientofthea1(j)termis30timeslessthan strong asymmetry between summer and winter on one
the observed difference between mean Dst values of side and spring and autumn on the other. To take it
summer and winter subsamples. This canbe explained into account we introduced additional terms into our
in the following way: there are other regressors, which regression, which are powers of a1(j) and their prod-
dependonparameterswithstatisticallysignificantsum- ucts with powers of a2(j). The sum of regressors with
mer/winterasymmetry,e.g. previousDstvalues. They the corresponding coefficients, depicted on Fig. 6, is
provide the lion share of summer/winter asymmetry of very similar to Fig. 4. Note that Fig. 6 was obtained
Dst. Agoodexampleofsucharegressoristhe sunspot independently from Fig. 4.
numberR,whichdescribesthe27-dayperiodicity. Nev- We did the same thing with the diurnal asymmetry.
ertheless, there is a small difference which can not be The distribution is plotted on Fig. 5, and the sum of
expressed with these terms. Including it into regres- regressorswith the corresponding coefficients – on Fig.
sion,weobtainthese statisticallysignificantregressors. 7. The term a1(j)·b1(j) is also significant and should
Tofurtherillustratethispoint,letusconsiderasanex- be included in the regression. After this we obtained a
ampleavalueX =const+Asinωt. Intheregressionit joint distribution of semiannual and diurnal variations
will look like Xn+1 =Xn+A[sinω(t+∆t)−sinωt]= of Dst index, plotted on Fig. 8. It contains 18 regres-
Xn+A[(cosω∆t−1)sinωt+cosωtsinω∆t]. Thefirst sors. Increasing the number of regressors describing
8
50
0.5
Prediction for
0 nts0.4 1 hour
poi 9 hours Bin size = 5 nT
Dst, nT-1-0500 on of all 00..23
cti
a
-150 Dst from Kyoto WDC for Geomagnetism Fr0.1
1h prediction
-200 1h prediction (autoregression) 0
-50 -40 -30 -20 -10 0 10 20 30 40 50
4-Nov 6-Nov 8-Nov 10-Nov12-Nov14-Nov16-Nov18-Nov20-Nov22-Nov Deviation, nT
1998 UT
Fig. 16 Error charts of predictions results of O’Brien &
Fig. 15 Comparison between predictions results of McPherron (2000a) (top) and ours (bottom) 1 hourahead.
O’Brien & McPherron (2000a) (top) and ours (bottom) 1 Errorchartforour9hourpredictionisplottedforreference
hour ahead. The following designations are used: ’Kyoto’
– official Dst index, available at Kyoto WDC for Geomag-
in much larger number of points with deviations more
netism;’AK1’–predictionbasedonthemodelofBurtonet
al.(1975)withre-calculatedcoefficients;’UCB’–prediction thanσ thanfollowsfromthe normaldistribution. This
based on the model of Fenrich & Luhmann (1998); ’AK2’ is mostly caused by the number of points in the wing
– prediction based on the model of O’Brien & McPherron bins |θ −hθ i| > 4σ being equal to 196 points ver-
V V
(2000b); ’ACE Gaps’ refer to the top line, indicating the sus 11 points in the case of normal distribution. How-
availability of solar wind data measured byACE satellite
ever, most of these points were obtained in the 1960s,
whenqualityofmeasurementswasmuchworsethento-
temporalvariationsofgeomagneticactivity we canim- day. This period includes the maximum and minimum
◦ ◦
prove the accuracy of this distribution. In particular, values of θV, equal to −59.7 and 18.8 . Neverthe-
one could add 11-year and 22-year solar cycles, higher less,these points constitute only aminor fractionofall
powers of a (j) and b (j) etc. points and didn’t affect the linear regression routine.
i i
Thus,we demonstratedhow easilyone cantake into Assumingnormaldistributionweobtainσ =2.925and
account new geoeffective parameters in this method’s hθVi =0.27<0.1σ. Thus, the distribution is insignifi-
framework. cantly shifted towards positive values.
Now let us discuss parameters, whose geoeffective- Ifweignorethewingbinsinthedistributionofmean
ness was determined by this method, and demonstrate DstvaluesagainstθV,whicharesomewhatrandomdue
that they are indeed geoeffective. to small amount of points in them, we will notice a
Latitudinalflowangleθ wasmostlyassociatedwith slight almost linear trend. If we plot the sum of terms
V
the southern component of IMF. I plotted the distri- containing θV (Fig. 10), we will notice a similar trend.
bution of its value and the corresponding mean Dst If we select two subsamples, one −8 < θV < −5 and
value on Fig. 9. The distribution looks similar to other 4 < θV < 9, and verify the hypothesis that the
a normal distribution, but it significantly differs from difference between the corresponding average Dst val-
the normal one according to χ2 test. This manifests uesis statistically significantusinga one-sidedStudent
Regressionmodelingmethod 9
50
0
-50
-100
T
n
Dst, -150 Dst from Kyoto WDC for Geomagnetism
-200
1h prediction
-250 1h prediction (autoregression)
-300
0
-350
-50
15-Jul 16-Jul 17-Jul 18-Jul 19-Jul 20-Jul -100
DOY 2000 -150
50 T-200
n
-500 Jul D1s5t- f1ro6m, 2K0y0o0to WDC Dst, --320500 Dst from WDC for Geomagnetism (Kyoto)
1h (autoregression) -350 1-hour forecast
T-100 -400 3-hour forecast
n
st, -150 -450
D-200 -500
-250 17-Nov 19-Nov 21-Nov 23-Nov 25-Nov 27-Nov
2003 UT
-300
-350 Fig. 18 Comparison between predictions results of Pal-
12:00 15:00 18:00 21:00 0:00 3:00 6:00 9:00 12:00 locchiaet al.(2006) (top)and ours(bottom)1hourahead.
UT hour
Our3-hourpredictionisgivenforreference. Onthetopplot
black line depicts Dst from Kyoto WDC, and the blue line
Fig.17 ComparisonbetweenpredictionsresultsofCerrato
– 1 hourprediction
et al. (2004), Burton et al. (1975), O’Brien & McPherron
(2000a),Fenrich&Luhmann(1998)(top)andours(middle)
1 hour ahead. Bottom plot is a scaled up version of the
Longitudinal flow angle ϕ was only occasionally
V
middle one
used in models. However, it appeared to be even more
significant than the latitudinal flow angle. Its distri-
test, we obtain t∞ = 6.278, which is over 99.95% sig- bution together with correspondingmean Dst values is
nificant. plotted on Fig. 12, where white bars show the com-
If we divide the sample in two subsamples: one for plete dataset sans rejects, and the grayed bars show
northern summer and one for northern winter, we ob- the quiet-time sample with Dst > −50nT. Like the
tain such a picture (Fig. 11). Once again, let us not latitudinalflowangle,the distributionofthe longitudi-
look at the wing bins. What do we see? Summer dis- nalflowangleresemblesnormaldistribution. However,
tribution has an obvious linear trend, but the winter χ2 test disproves the relevant null-hypothesis. Once
one has not. If we apply the Student test to the same again,this is mostlydue to wing bins whicharemostly
intervals now, we obtainthat t∞ is 5.44 in the summer formed of data points, corresponding to measurements
and only 0.059 in the winter. The prior corresponds to in 1960s, including maximum and minimum values
morethan99.95%significance,whilethelatter–toless equal to −65.6◦ and 48.5◦. Assuming normal distribu-
than 10%. This could mean that there are two factors tion we obtain σ =2.934 and hϕ i=−0.30≈−0.1σ.
V
connected with the latitudinal flow angle, which work A significant trend is the most prominent feature of
together in the summer and against each other in the this figure. If we plot a sum of regressors, which con-
winter. The physical explanation of this phenomenon, tain ϕ (Fig. 13), we see a very similar trend. Like
V
however, lies beyond the scope of this paper. before, we plot the distribution for summer and winter
10
subsamplesseparately(Fig. 14). Weseethatthetrend (“test” sample) or not. To further illustrate this point
is identical on both plots, so the corresponding effect let us consider an example: I predicted Dst 3 hours
is season-independent. The list of regressors for k =1, ahead using the “test” sample from 2007 to 2008 and
containing θ and ϕ , is given in Table 1. It contains 3 “training” samples: the sample from 1976 to 2003
V V
the regressorsthemself, their coefficients and F values. gives LC=0.851, the sample from 1976 to 2006 gives
Thus,we demonstratedthat ourmethod is truly ca- LC=0.854, and the sample from 1976 to 2008 gives
pable of pointing out new geoeffective parameters and LC=0.854 as well. The first two “training” samples
verified the geoffectiveness of two such quantities. do not contain the “test” sample, but the results are
slightlydifferent. The third“training”samplecontains
the “test”sample, yet the result is the same as for the
6 Prediction results second sample. This yields a conclusion that the vol-
ume and statisticalproperties of the “training”sample
Takingintoaccounttheconsideredparameterstogether affectthe predictionresultsstrongerthanthe inclusion
with parameters, whose geoeffectiveness was beyond of the “test” sample. So, the inclusion of the “test”
doubts, like previous values of Dst, dawn-dusk electric sample to the “training” sample has little or no effect
field, ram pressure of the solar wind and most of other on the LC. PE is calculated independently from the
parameters from the OMNI2 database, we constructed “test”sampleandisnotaffectedbyits selectioninany
models for predicting Dst 1, 3, 6, 9, 12, 18, and 24 way.
hours ahead, and 3 more models for predicting Dst 1 We also present graphical representations of pre-
hour ahead for quiet and perturbed conditions and for diction in comparison with results of other authors:
thecasewhensatellitedataareunavailable(autoregres- O’Brien & McPherron (2000a) (Fig. 15, 16), Cid et
sion, see more in Parnowski (2009b)). The statistical al. (2005) (Fig. 17), and Pallocchia et al. (2006) (Fig.
characteristicsofthesemodelsaresummarizedinTable 18). Note,thattheintervalsforpredictionwereselected
2. They include Residual Mean Square (RMS), Linear by authors of the original papers, and the coefficients
CorrelationCoefficient (LC), and Prediction Efficiency in eq. (1) were the same all the time and for all the
(PE =1−RMS2/SD2,whereSDisthesample’sStan- figures. It is clearly visible that our method provides
dard Deviation). In divided cells the top number cor- much more precise forecast than most empirical mod-
responds to the actual model, and the bottom one – els and typical neural network models. Ridiculously,
to the simplest possible model D (j+k)=D (j). It even our 9h model is more precise than most empir-
st st
is noteworthy that despites good correlationfor all the ical 1h models. The autoregression model, described
models, in reality only the 1-hour and 3-hour models by eq. (2), though, lags in the left part of the plot
are ready for practical use, and the 6-hour model can due to a rapid positive change of Dst at 1500 UT. The
potentially reachthis state. This is due to a significant lag persists through the growth phase and the main
time shift being present in further predicting models. phase, and vanishes only in the recovery phase. For
Note that since the proposed method is statisti- thisreason,theautocorrelationmodelholdslittleprac-
cal, there is little difference whether the “training” tical value and should be considered as a transitional
sample contains the period when prediction is made result, required to construct the full model. It is, how-
ever,possible to improve it by adding terms describing
temporal variations, and, for example, the number of
Table1 Regressorscontainingtheflowangles. V(j)isthe sunspots, but then the term “autoregression” will no
bulk flow velocity of thesolar wind longer be applicable.
i x C F On Fig. 19 we present the results of prediction 3,
i i i
1 θ (j)·V(j) (2.8±0.9)·10−5 10.2 6 and 9 hours ahead for a number of events, kindly
V
2 θ4(j) (−4.9±1.2)·10−5 17.5 selected for us by V.G. Fainshtein, which are particu-
V
3 a1(j)·θV(j)·Dst(j) (1.3±0.3)·10−3 18.7 larly hard to predict by medium-term methods, such
4 a1(j)·θV4(j) (−3.2±1.6)·10−5 4.0 asEselevichet al.(2009),to verify the efficiency of our
5 ϕ (j) (−4.1±0.6)·10−2 50.1 method. We can see that this method’s accuracy is
V
6 ϕ2(j) (−4.5±1.5)·10−3 9.3 higher for stronger storms, which are of greater inter-
V
7 ϕ3(j) (−3.5±1.3)·10−4 7.1 est. A huge advantage of this method is that the most
V
8 ϕ4(j) (5.3±1.7)·10−5 10.1 resource-demanding operation – the calculation of the
V
9 a2(j)·ϕV(j) (−2.2±0.8)·10−2 7.9 regressioncoefficients – should be performed only once
10 a2(j)·ϕ3V(j) (5.1±1.9)·10−4 7.6 for each model. The prediction itself is just a summa-
11 a2(j)·ϕ4V(j) (−3.5±1.4)·10−5 6.6 tionofapolynom,whichusuallytakesnomorethan4-6