# Quasi-Periodic Peak Energy Oscillations in X-ray Bursts from SGR J1935+2154

Oliver J. Roberts Science and Technology Institute, Universities Space and
Research Association, 320 Sparkman Drive, Huntsville, AL 35805, USA.  Matthew
G. Baring Department of Physics and Astronomy - MS 108, Rice University, 6100
Main Street, Houston, Texas 77251-1892, USA  Daniela Huppenkothen SRON
Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The
Netherlands  Ersin Göğüş Sabancı University, Faculty of Engineering and
Natural Sciences, İstanbul 34956 Turkey  Yuki Kaneko Sabancı University,
Faculty of Engineering and Natural Sciences, İstanbul 34956 Turkey  Chryssa
Kouveliotou Department of Physics, The George Washington University, 725 21st
Street NW, Washington, DC 20052, USA  Lin Lin Department of Astronomy, Beijing
Normal University, Beijing 100875, China  Alexander J. van der Horst
Department of Physics, The George Washington University, 725 21st Street NW,
Washington, DC 20052, USA  George Younes NASA Goddard Space Flight Center,
University of Maryland, Baltimore County, Greenbelt, MD 20771, USA

###### Abstract

Magnetars are young neutron stars powered by the strongest magnetic fields in
the Universe (1013-15 G). Their transient X-ray emission usually manifests as
short (a few hundred milliseconds), bright, energetic ($\sim$ 1040-41 erg)
X-ray bursts. Since its discovery in 2014, magnetar SGR J1935+2154 has become
one of the most prolific magnetars, exhibiting very active bursting episodes,
and other fascinating events such as pulse timing anti-glitches and Fast Radio
Bursts. Here, we present evidence for possible 42 Hz (24 ms) quasi-periodic
oscillations in the $\nu F_{\nu}$ spectrum peak energy ($E_{\rm p}$)
identified in a unique burst detected with the Fermi Gamma-ray Burst Monitor
in January 2022. While quasi-periodic oscillations have been previously
reported in the intensity of magnetar burst lightcurves, quasi-periodic
oscillations in the $E_{\rm p}$ have not. We also find an additional event
from the same outburst that appears to exhibit similar character in $E_{\rm
p}$, albeit of lower statistical quality. For these two exceptional
transients, such $E_{\rm p}$ oscillations can be explained by magnetospheric
density and pressure perturbations. For burst-emitting plasma consisting
purely of $e^{+}e^{-}$ pairs, these acoustic modes propagate along a highly
magnetized flux tube of length up to around $L\sim 130$ neutron star radii,
with $L$ being lower if ions are present in the emission zone. Detailed time-
resolved analyses of other magnetar bursts are encouraged to evaluate the
rarity of these events and their underlying mechanisms.

Neutron stars (1108); Magnetars (992); Compact objects (288); Soft gamma-ray
repeaters (1471)

##  1 Introduction

SGR J1935+2154 (hereafter SGR J1935) was discovered in 2014 (Stamatikos et
al., 2014) with the Burst Alert Telescope (BAT) on the Neil Gehrels Swift
Observatory (hereafter Swift), when the source emitted a few short X-ray
bursts. Further observations of the source with the Swift X-ray Telescope
(XRT), Chandra and XMM-Newton observatories, found a spin period of 3.24 s and
a period derivative of 1.43$\times$10-11 s$\cdot$s-1, indicating a dipole
magnetic field strength of $\sim$2.2 $\times$ 1014 G (Israel et al., 2016) at
the equator, which confirmed the magnetar nature of the source. The source has
exhibited prolific bursting behavior since its discovery, with multiple burst-
active episodes in 2015 and 2016 (Younes et al., 2017), 2018, 2019 and 2020
(Lin et al., 2020), and in 2021 and 2022 (Roberts, 2023). Comprehensive
investigations of SGR J1935 bursts, using data from the $Fermi$ Gamma-ray
Space Telescope’s Gamma-ray Burst Monitor (Fermi/GBM) and Swift/BAT during
these active episodes, showed that the magnetar became progressively more
burst-active in every subsequent outburst episode (Lin et al., 2020). The 2020
SGR J1935 outburst was the most energetic, producing a “burst storm” on April
27th detected with multiple instruments, including Fermi/GBM (Kaneko et al.,
2021) and NICER (Younes et al., 2020), the latter reporting a burst rate
maximum of 0.2 bursts/s over a $\sim$1200 s pointed observation. A hard burst
with a $\nu$Fν peak energy ($E_{\rm p}$) reaching up to $\sim$80 keV on April
28th, which was contemporaneous with a Fast Radio Burst (FRB; Ridnaia et al.,
2021; CHIME/FRB Collaboration et al., 2020; Bochenek et al., 2020; Mereghetti
et al., 2020), confirmed that at least one magnetar, SGR J1935, was a source
of these exotic radio transients. SGR J1935 also underwent a radio-pulsar
phase in October 2022, a few days after an anti-glitch (i.e., sudden spin
down) was inferred from X-ray monitoring observations (Younes et al., 2022).
Two additional FRB-like bursts from the source were detected during the same
period (Kirsten et al., 2021).

In this paper, we present $E_{\rm p}$ oscillations found in two bursts from
SGR J1935, detected on the 12th and 16th of January, 2022. These observations
open up new territory for probing the physics of magnetar bursts and the
character of neutron star magnetospheres. In Section 2, we introduce the data
analysis methods and spectral analysis. In Section 3, we identify through
visual inspection of the brightest bursts from SGR J1935 ($f$ $\geq$
1$\times$10-6), two candidates that appear to exhibit $E_{\rm p}$
oscillations. We present the time-integrated spectral fits for each burst, the
time-resolved energy spectra, a periodicity search and temporal analysis. We
present a theoretical interpretation of these results in Section 4, and
summarize our findings in Section 5.

##  2 Data Analyses

We use the triggered Time-Tagged Event (TTE) data (time resolution of $\sim$2
$\mu$s, 128 pseudo-logarithmically spaced energy channels), obtained with the
Fermi Gamma-ray Burst Monitor’s (GBM’s) thallium-doped Sodium Iodide (NaI)
detectors with good detector-to-source angles ($\leq$60∘). These detectors
have an effective spectral range of $\sim$8$-$900 keV, which nicely overlaps
with the energy range of normal magnetar burst spectra ($\leq$300 keV). More
information on the Fermi/GBM instrument and its data types can be found in
Meegan et al. (2009) and von Kienlin et al. (2020).

###  2.1 Spectral Analysis

We performed spectral analysis on bright SGR J1935 bursts, emitted during the
source activation starting in 2019 through 2022. While analyzing the Fermi/GBM
TTE data, we found two bursts (labeled throughout this study as “Burst 1” and
“Burst 2”), that appear to show variations in their $\nu$Fν peak or $E_{\rm
p}$ values. We performed a time-resolved spectral analysis of these bursts
using a power-law with an exponential cutoff, commonly referred to as a
Comptonized (COMPT) model. This model is defined as:

$F(E)=A\Biggl{(}\frac{E}{E_{\rm
piv}}\Biggl{)}^{\alpha}\exp\left[-\frac{(\alpha+2)E}{E_{\rm p}}\right],$ (1)

where, _A_ is the amplitude in ph/s/cm2/keV, $E_{\rm p}$ is the $\nu$Fν peak
in keV, _$\alpha$_ is the power-law index, _$E_{\rm piv}$_ is the pivot energy
in keV, fixed at 100 keV in this study.

In addition to performing time-integrated and time-resolved spectral analysis
using the COMPT model, other models such as a simple power-law (PL), an
Optically-Thin Thermal Bremsstrahlung (OTTB), a single black-body (BB) and a
double black-body (BB+BB), were also used. These were fit to the burst spectra
using RMFIT (v4.4.2), developed specifically for the analysis of GBM
data111https://fermi.gsfc.nasa.gov/ssc/data/analysis/rmfit. We selected
background intervals devoid of bursts several seconds before and after each
event, fitting each interval with a polynomial function to model the
background during the burst, which was subsequently subtracted. The Detector
Response Matrices (DRMs) for triggered bursts are available as a part of the
publicly accessible data products. For untriggered bursts, we generated DRMs
using GBMRSP v2.0. For the spectral analysis of all the bursts in our sample,
we used the NaI-detector energy range of $8-900$ keV. We neglect an energy
interval of $\sim$4 keV centered around 35 keV, where the k-edge feature from
Iodine in the NaI detectors appears in the data, as this has not been modeled
perfectly (Bissaldi et al., 2009). The exclusion of this small portion of data
from our analysis does not affect the results, but improves the statistics.

We use the Bayesian Information Criterion (BIC; Schwarz, 1978; Liddle, 2007)
to determine preferred models for the bursts in this study, which is often
used for model comparisons with the maximum likelihood statistics. We
calculated BIC for each spectral fit as follows:

${\rm BIC}=-2\ln\mathcal{L}_{\rm max}+k\ln N={\rm CSTAT}+k\ln N,$ (2)

where $\mathcal{L}_{\rm max}$ is the maximum likelihood, $k$ is the number of
free parameters in the spectral model and $N$ is the number of data points. We
then calculated $\Delta$BIC for each pair of the four models for comparing the
posterior probabilities of the given two models. We employ $|\Delta$BIC$|\geq
12$ to select a preferred model, which corresponds to the Bayes factor of
$\sim$400, indicating that the posterior probability of the model with a
smaller BIC value is higher by $>99.7\%$ (Kass & Raftery, 1995; Anderson &
Burnham, 2004; Liddle, 2007).

We only report the parameters from the COMPT and BB+BB model fits to our data,
as they fit the data better than the PL, OTTB and BB models. We compare the
BIC values between the COMPT and BB+BB model fits by defining $\Delta$BIC =
BICCOMPT – BICBB+BB, where a BB+BB model is the preferred fit to the burst
spectrum when $\Delta$BIC $\geq$ 12, and a COMPT model is the preferred fit if
$\Delta$BIC $\leq-12$, at 3-$\sigma$ significance. The duration and spectral
parameters from fitting the candidates using the aforementioned models, are
shown in Table 1.

###  2.2 Localization and Burst History

In the following sections we discuss the time-history and localization of the
two burst outliers detected in January 2022. The most recent event (Burst 1;
see Table 1) was localized to a Right Ascension (RA) and Declination
($\delta$) of 294.2∘ and 23.8∘ respectively, within the 1-$\sigma$ error
(3.6∘) of the known position of SGR J1935+2154 (RA: 293.7∘, $\delta$: 21.9∘
(Israel et al., 2016)). Similarly, the other event (Burst 2, which occurred
around 3.8 days earlier than Burst 1) was localized to RA: 291.8∘, $\delta$:
24.5∘, with an error of 4.5∘ (1-$\sigma$), 3.3∘ degrees from the SGR
J1935+2154 position. Both events are consistent, within their errors, with an
origin from SGR J1935; they occurred during a heightened period of burst-
activity from the source, detected by multiple instruments (i.e., Roberts &
Fermi GBM Team (2022); Ridnaia et al. (2022); Kozyrev et al. (2022)). We used
a Bayesian block algorithm, similar to that used in Lin et al. (2020), to
determine the duration of each burst ($\tau_{BB}$), shown in Table 1.

##  3 Results

The time-integrated COMPT and BB+BB spectral fits for Burst 1 and Burst 2 are
presented in Table 1, noting that the spectral analysis for Burst 1 was
performed separately for the whole event (which includes the faint, weak
emission), as well as solely for the brighter, second emission episode. The
fits were performed over the energy range $\,10-1000\,$keV, though most burst
counts above background were at energies $\,\lesssim 100\,$keV. We find that a
COMPT model fit to the spectrum is preferred for both candidates.

###  3.1 SGR J1935 $E_{\rm p}$ Oscillation Candidates

The best example of fluctuations in $E_{\rm p}$ is Burst 1, which triggered
Fermi/GBM at 14:09:38.94 UTC on January 16th, 2022, and subsequently assigned
the burst number, bn220116590. The burst appears in the data as a faint pulse
followed about 100 ms later by a very bright pulse, shown in the left panel of
Fig. 1, with a total duration of 693 ms. The time-integrated spectral results
for both the whole burst and the brighter, second emission episode are shown
separately in Table 1. We use NaI detectors 0, 1, 2, 3 and 5, as these all
have good detector-source angles ($\leq$60∘). Upon fitting the brighter,
second pulse, we find that the overall $E_{\rm p}$ has a soft-hard-soft trend
over that entire emission duration. When using finer temporal binning, the
spectrum suggests possible periodic perturbations in $E_{\rm p}$ (right panel
of Fig. 1). Using the distance to the supernova remnant G57.2+0.8 (purported
to be the host of SGR J1935), of 9 kpc (Pavlović et al., 2013; Zhong et al.,
2020), we derive Eiso and Liso values of 3.9$\times$1040 erg and
1.9$\times$1041 erg/s, respectively for the second, 200 ms emission episode.

The next best example of an event with fluctuations in $E_{\rm p}$ is Burst 2,
which triggered Fermi/GBM at 19:58:04.04 UTC on January 12th, 2022, and was
subsequently assigned the burst number, bn22012832. Burst 2 appears as a
bright burst with a duration of 254 ms. We use NaI detectors 0, 1, 3 and 5, as
these all had good detector-source angles ($<$60∘). When fitting a COMPT model
to the spectrum, the general $E_{\rm p}$ behavior of the burst is a quick rise
followed by a slow decay, with weak perturbations on a similar timescale to
that observed in Burst 1 (see Fig. 1). The time-integrated spectral results
for this burst are shown in Table 1. Using the distance to SGR J1935 of 9 kpc,
we derive Eiso and Liso values of 2.2$\times$1040 erg and 8.8$\times$1040
erg/s, respectively for this 254 ms event.

Table 1: Ep Oscillation Candidates

ID | Date | Trigger Time | $\tau_{BB}$ | Fitted Interval (ms) | Comptonized (COMPT) | BB+BB |  $\Delta$BIC | Fluence
---|---|---|---|---|---|---|---|---
| YYMMDD | UTC | (ms) | tstart:tstop | $\alpha$ |  $E_{\rm p}$ (keV) |  $kT_{1}$ (keV) |  $kT_{2}$ (keV) |  | ($\times$10-6 erg cm-2 )
1 | 220116 (i) | 14:09:38.94 | 693 |  $-$88 : 320 | 0.38$\pm$0.05 | 33.0$\pm$0.2 | 5.5$\pm$0.2 | 10.8$\pm$0.2 |  $-$16 | 1.95$\pm$0.01
| 220116 (ii) | - | - | 120 : 320 | 0.47$\pm$0.05 | 34.0$\pm$0.2 | 5.9$\pm$0.2 | 11.0$\pm$0.2 | - | 3.98$\pm$0.03
2 | 220112 | 19:58:04.04 | 254 |  $-$24 : 276 | 0.51$\pm$0.07 | 30.1$\pm$0.3 | 4.7$\pm$0.3 | 9.3$\pm$0.2 |  $-$11 | 2.30$\pm$0.02

Note. — The time-integrated spectral analysis for Burst 1 was performed both
over the whole burst (i) and for the second, brighter emission (ii). All
statistical uncertainties are at the 1$\sigma$ confidence level. The tstart
and tstop are relative to the trigger time (in UTC). The fluence values are
over $\,10-1000\,$ keV using a COMPT spectrum.

Comparing both burst spectra, the burst spectrum of Burst 1 is found to have a
clear ingress and egress in $E_{\rm p}$ overall (i.e., the envelope), peaking
at around 50 ms after the onset of the second emission episode (25 ms in the
right panel of Fig. 1)222This soft-hard-soft evolution is suggestive of a
radiation cone passing into and out of the line-of-sight of Fermi/GBM.. This
is markedly different from the flatter overall $E_{\rm p}$ behavior of Burst 2
and likely implies different viewing geometries of the source environment.
While a BB+BB model was not found to be the preferred fit to either burst
spectrum, we calculate the BB region sizes to try to put constraints on the
flux tube geometry. For Burst 1, we find the BB regions to be 283 $\pm$ 52 km2
and 31 $\pm$ 3 km2 for $kT_{1}$ and $kT_{2}$, respectively. Similarly, we
calculate the BB regions for Burst 2 to be 359 $\pm$ 125 km2 and 57 $\pm$ 8
km2 for $kT_{1}$ and $kT_{2}$, respectively. The $kT_{1}$ region size is used
to provide a constraint on the diameter of the cross section of the activated
flux tube, as described in Section 4.3.

The light curve and $E_{\rm p}$ time traces in Figure 1 display similar
character in their general shapes, or envelopes, suggesting that there is a
significant correlation between the flux $\cal{F}$ and $E_{\rm p}$. To explore
this, we binned flux and $E_{\rm p}$ data in 4 ms time intervals spanning the
burst from the rapid risetime and decay of each event. A correlation of ${\cal
F}\,\propto\,E_{\rm p}^{3.2\pm 0.2}$ was found, with a Spearman rank order
correlation coefficient and chance probability of $0.91$ and $5.50\times
10^{-29}$, respectively. This differed considerably from the ${\cal
F}\,\propto\,E_{\rm p}^{2}$ correlation obtained by (Roberts et al., 2021) for
the GRB 200415A magnetar giant flare from the Sculptor galaxy, a correlation
that is nominally a signature of Doppler beaming/boosting from ultra-
relativistic outflows from a rotating star. Thus, we conclude here that the
plasma that generated these two special bursts was likely moving only mildly-
relativistically along the active field lines, with the radiation being only
somewhat anisotropic at each emission locale.

Figure 1: Left: The lightcurves of Burst 1 (bn220116590) and Burst 2
(bn22012832) plotted alongside each other within a 500 ms window starting at
-130 ms. The trigger time for each event is the zero time. The first, weaker
emission episode in Burst 1 lasts about 100 ms, starting at $\sim$125 ms.
Right: A $\sim$175 ms window showing the $E_{\rm p}$ behavior of Burst 1
(pink) and Burst 2 (purple). Burst 1 is shifted in time by -166 ms in order to
show the similar $E_{\rm p}$ oscillating frequency of 24 ms (42 Hz),
highlighted by the dashed grey lines. The temporal binning for both panels is
4 ms.

###  3.2 Periodicity Search for Burst 1

We performed three complementary searches of the time series of $E_{\rm p}$
values of Burst 1 at both $2~{}\mathrm{ms}$ and $4~{}\mathrm{ms}$ resolution
for statistical evidence of a periodic or quasi-periodic oscillation (QPO);
two Fourier approaches and a more sophisticated Gaussian Processes technique.
As Burst 2 was found to be weaker than Burst 1, we did not search its time
series for any spectral QPOs (abbreviated hereafter QPSOs333In order to
differentiate between traditional QPOs and the spectral QPOs found in this
study, we named the latter with the abbreviation, QPSO (quasi-periodic
spectral oscillations) throughout the remainder of the paper.) in this event.
To search for periodic signals in the $E_{\rm p}$ lightcurves, we follow the
procedure of Huppenkothen et al. (2013). We first generate a periodogram from
the $E_{\rm p}$ time series, fit a power law period distribution to this
periodogram, and subsequently find the highest maximum outlier. The
significance is calibrated using $1000$ simulated periodograms based on
posterior distributions of the power-law model parameters generated using
Markov Chain Monte Carlo simulations. We assumed wide, flat priors for the
power-law index and the logarithm of the power-law amplitude, and a normal
prior with a mean of $2$ and a width of $0.5$ around a noise value of $2$. We
find no statistical evidence to reject the null hypothesis, namely that the
periodogram is consistent with stochastic variability following a power-law
power spectrum, and thus find no evidence for a periodic signal in the data
using this analysis.

Next, we search for QPSOs by comparing the power-law model for the periodogram
to a model consisting of a power-law model and a QPSO – parametrized as a
Lorentzian – through a Likelihood Ratio Test (LRT). Once again, the LRT is
calibrated via $1000$ simulations generated from an MCMC sample of the null
hypothesis (the power-law model). Here, we find weak evidence in the
$2\,\mathrm{ms}$ data ($p=0.039$) that the observations might not be drawn
from the null hypothesis. As the $E_{\rm p}$ values are spectral model
parameters associated with uncertainties, and because Fourier periodograms by
default do not take uncertainties into account, we check the effects of
incorporating the uncertainties through simulations. Here, for each point in
the time series, we draw from a normal distribution with the mean given by the
best-fit $E_{\rm p}$ value, and the standard deviation given by the $1\sigma$
parameter uncertainties. In this way, we draw $1000$ time series.
Subsequently, we compute the LRT between the models with and without a QPSO
for each simulated time series, and compare the resulting distribution to the
distribution generated from the MCMC simulations of only the power-law model.
We find that, overall, while not exactly the same, the LRT of simulations
generated from the best-fit values and their uncertainties have a similar
distribution to the LRTs generated from the power-law model only, with no QPSO
present. The periodogram of the best-fit values is a mild outlier ($p\sim
0.04$) with respect to both distributions, and thus we conclude that there is
weak evidence for a rejection of the null hypothesis that the periodogram was
generated by pure noise.

Figure 2: Searching for (quasi-)periodicities in $E_{\rm p}$ with a Gaussian
Process-based model. Left: Ep values with uncertainties (black) as a function
of time, along with posterior draws of the skew-Gaussian mean function (green)
and posterior draws of the combined model of skew-Gaussian and damped random
walk in orange. Overall, this model, containing only aperiodic stochastic
variability on top of the overall skew-Gaussian trend, fits the $E_{\rm p}$
very well. Middle: as in the left panel, but for a model where the damped
random walk was replaced by a QPSO, parametrized as a stochastically driven,
damped harmonic oscillator. This model, too, describes the data adequately
well. The resulting Bayes factor does not allow us to favour one model over
the other. Right: the posterior probability density for the period of the QPSO
in the model shown in the middle panel, showing that the time series
constrains the QPSO period well.

Fourier-based methods have a number of shortcomings in this context. Firstly,
they assume a stationary process, which the $E_{\rm p}$ time series of a
magnetar burst is not, and which is known to bias the significance of QPO
detections (Hübner et al., 2022a). Secondly, as explained above, they cannot
robustly take into account the uncertainties on the $E_{\rm p}$ value.
Finally, null hypothesis tests can only ever reject a null hypothesis, and
thus cannot yield direct evidence for the detection of a QPSO. With this in
mind, we implement the Gaussian Process-based method recently introduced in
Hübner et al. (2022b). Gaussian Processes encompass a class of time domain-
based models that addresses many of the above shortcomings. They enable direct
(Bayesian) model comparison via comparison of model evidences in the form of
Bayes factors. They can natively include non-stationarity via mean functions,
and they include prescriptions to take into account data uncertainties. We
choose a skew-Gaussian function as a mean model to account for the broad trend
in the time series, and implement three different hypotheses for the
variability seen on shorter timescales than that trend: a damped random walk
in the $E_{\rm p}$ variate to parametrize stochastic, aperiodic variability; a
stochastically driven, damped harmonic oscillator as a QPSO model
(corresponding to a Lorentzian in the Fourier domain); and a model combining
both damped random walk and Lorentzian.

For all models, we implement wide, uninformative priors as described in Hübner
et al. (2022b), sample the posterior of each model via Dynamic Nested Sampling
(Speagle, 2020), and finally compute Bayes factors to compare all models.
Unlike null hypothesis tests, Bayes factors directly compare two models to
each other, and thus can be used to yield evidence for or against a given
model. We find that the QPSO-only and the red noise-only models produce very
similar evidences, $\log(B)=0.018$ for both $2~{}\mathrm{ms}$ and
$4~{}\mathrm{ms}$ time series. This is smaller than the uncertainties in the
calculation of the evidences $\sigma_{B}\sim 0.2$, and both models thus
describe the data equally well. We note that both are mildly preferred over
the model containing both red noise and a QPSO, $\log(B)=2$, likely owing to
the increase in model parameters (and thus prior volume) of the more complex
model. In Figure 2, we present the $2~{}\mathrm{ms}$ time series along with
posterior draws from the model containing only the mean function and red
noise, and for the model containing only the mean function and a QPSO. We find
that indeed, both are able to describe the observed data well. We also note,
however, that the posterior probability density for the QPSO’s period
parameter is narrow and unimodal, $P=0.02415^{+0.0023}_{-0.0018}$ seconds,
suggesting that the QPSO is well-constrained.

##  4 Discussion

To recap, the time-integrated spectra of the two selected bright bursts from
SGR J1935+2154  were found to prefer a non-thermal fit, and values of $E_{\rm
iso}$ and $L_{\rm iso}$ for these two transients were comparable to those of
typical energetic magnetar bursts. Our time-resolved analysis of $E_{\rm p}$
in the brighter event (Burst 1) uncovered suggestive evidence of $E_{\rm p}$
fluctuations, contrasting its lightcurve morphology which doesn’t exhibit
traditional QPO signatures. Careful and thorough analysis to evaluate the
periodicity in these $E_{\rm p}$ oscillations uncovered a well-constrained, 24
ms (i.e., 42 Hz) QPSO, the first time such possible periodic spectral
fluctuations have been observed. Here we address some nuances of the QPSO
analysis, summarize prior evidence for light curve QPOs in bursts, and deliver
a theoretical interpretation for what could generate QPSOs in these two
events.

###  4.1 Bayes Factors in QPSO Probes

It is important to mention that there are several caveats concerning the Bayes
factors adopted in the $E_{\rm p}$ oscillation analysis. First, the models
that were compared were purely empirical choices, and in the absence of a
physical model capable of generating the data we see, were our only option for
searching for QPSOs in this kind of data. The damped random walk model makes
strong assumptions about the data-generating process; it is worth bearing in
mind that its standard assumptions originate in studies of black holes, not
magnetars. We expect that a more physically-motivated model for generating the
variability would render the Bayes factor a more reliable metric for model
selection. Second, Bayes factors are known to be very sensitive to prior
choices. Here we set wide, uninformative priors on the model parameters, so
that our lack of physical knowledge of the data-generating process likely
impacts the significance of the result. Third, a statistically sound analysis
should rely on more than a single number. Here, we jointly considered the
evidence from both the Fourier-based analysis, which yielded a mild outlier,
and the Gaussian Process analysis, which yielded an inconclusive Bayes factor,
but a large effect size in the form of a narrow, well-constrained posterior on
the QPSO period. For a classical noise process, we would expect the period
posterior to largely reproduce the wide, flat, uninformative prior rather than
be heavily concentrated and unimodal, as is seen here.

Finally, it is important to note the impact of chosen priors on the models to
be compared, which in the Gaussian approach were the damped random walk and
the Lorentzian QPSO. In the analysis above, we have considered the Bayes
factor in isolation, which in practice is equivalent to assigning equal
probability to both the damped random walk and the Lorentzian model. If QPSOs
are physically expected, as discussed at length in Section 4.3, then priors
change and the confidence in the existence of the 24 ms QPSO in Burst 1 rises
accordingly. The combination of these considerations leads us to conclude that
this candidate discovery is worth reporting, especially in the context of the
observation and analysis of future bright bursts, where improved count
statistics and longer durations may yield more confident detections.

###  4.2 QPOs in Magnetar Bursts

Quasi-periodic oscillations with frequencies of $\sim$18-625 Hz have been
identified in light curves from Galactic magnetar giant flares from SGR
1806-20 and SGR 1900+14 (Strohmayer & Watts, 2005, 2006; Huppenkothen et al.,
2014b; Miller et al., 2019). Recently, QPOs have also been identified in the
tail (7 to 160 ms) of extragalactic magnetar giant flare GRB 200415A
(183$\pm$20 Hz at 2.5$\sigma$ signficance, p-value of
2.3$\times$10${}^{-2}\,$; see Roberts et al., 2021), and during the peak (3
ms) of GRB 200415A (2132 and 4250 Hz; see Castro-Tirado et al., 2021).
Huppenkothen et al. (2014c) identified a weak QPO signal centered at $\sim$260
Hz in an SGR J1550-5418 burst, and more significant QPOs centered at $\sim$93
and $\sim$127 Hz when using stacking of the same data. A QPO centered at
$\sim$57 Hz was also found by stacking 30 short, individual bursts from SGR
1806-20 (Huppenkothen et al., 2014a). More recently, a 40 Hz QPO at
3.4$\sigma$ (p-value of 2.9$\times$10-4) was reported from SGR J1935+2154 in
data from the Hard X-ray Modulation Telescope (HXMT; an instrument on Insight)
over an energy interval of 18-50 keV (Li et al., 2022). This unusual burst was
observed contemporaneously with FRB 200428 during the 2020 outburst,
suggesting that some FRBs are related to strong oscillation processes in
neutron stars.

All these QPOs were identified in the intensity of the lightcurve over
frequencies predominately from 10–2000 Hz. These are most likely associated
with global torsional/axial oscillations of the neutron star, particularly in
the lower frequency band below 150 Hz (Levin, 2007; Huppenkothen et al.,
2014b; Miller et al., 2019). However, models to explain spectral fluctuations
(QPSOs) at similar frequencies are lacking, likely primarily because such a
phenomenon has not been reported before. Consequently, we now present an
original theoretical interpretation of the QPSOs ($E_{\rm p}$ oscillations)
suggested by the spectro-temporal analysis above.

###  4.3 Magnetar Burst Acoustics

The emission region for the highly luminous bursts necessarily is highly
optically thick to Thomson scattering due to the high densities of the
radiating charges (e.g., Lin et al., 2011, 2012). Thus the microphysical
fluctuation (Thomson diffusion) timescales should be very small. In addition,
the Alfvén mode frequencies lie below the ion plasma frequency
$\omega_{p}=\sqrt{4\pi e^{2}n_{p}/m_{p}}$, which is likely $\sim 10^{14}\,$Hz
for plasma number densities of $n_{p}\sim 10^{22}\,$cm-3 that can be deduced
from burst radiation efficiency arguments (Baring & Harding, 2007). While very
low frequency (long wavelength) Alfvén waves can be contemplated, the high
plasma density effectively scatters and absorbs them. The result is a
decoherence that indicates that Alfvén modes do not drive the $E_{\rm p}$
oscillations that are clearly a macro-scale phenomenon. The high opacity
implies that the emission environment should approximate, to a considerable
extent, a quasi-thermodynamic ensemble of photons, pairs and possibly also
ions.

The natural manifestation of thermodynamic fluctuations in such a highly
magnetic environment is via density ($\,\rho\,$) perturbations along surface-
to-surface magnetic flux tubes accompanied by pressure ($\,P\,$) variations
mediated by the gaseous equation of state (EOS), i.e. acoustic oscillations.
Pressure changes/fluctuations will then generate varying values for the
observed photon $E_{\rm p}$ value. Generally, the anticipated chaotic nature
of the burst region should preclude clear oscillation signatures if the
emission zone straddles many field lines. Yet, occasionally, the activated
flux tube could possess a range of small cross sectional areas $A$, and then
pressure wave signatures might not be totally obscured by the environmental
chaos. The radiating plasma “bounces” between flux tube footpoints, with
adiabatic expansion and compression along the tube seeding the $E_{\rm p}$
oscillations. We explore here the implications of this picture.

A fairly “constrained” magnetic flux tube with a small cross section can serve
as an acoustic cavity for the burst zone, so that natural oscillations can be
permitted along its length on timescales of arclength $\,\cal S\,$ of the tube
divided by the sound speed $\,c_{s}\,$. This should be commensurate with the
favored $\,24\,$ms period in the QPO analysis of the $E_{\rm p}$ variability.
Observe, that the discussion above indicated that the historical explanations
of QPOs in the light curves of magnetar giant flares (e.g., Strohmayer &
Watts, 2005, 2006) have considered sub-surface seismic modes. In principle,
these could provide a driver for the $E_{\rm p}$ fluctuations. Yet, any
periodic imprint of sub-surface drivers for burst activation is generally
obscured by the high opacity of the flux tube, unless the seismic mode
frequency approximately coincides with the natural frequency of the acoustic
cavity, in which case the sub-surface driving is resonant. This may or may not
be the case for the bursts studied here.

Returning to the acoustics, since the radiation efficiency of the burst region
is likely small, and the photon peak energies are well below the electron rest
mass energy, the system pressure is putatively dominated by the gas
contribution. The sound speed can then be estimated using the formalism
provided by Synge (1957) for a 3D relativistic Maxwell-Boltzmann distribution
of a single species. For the purposes of this discussion, we will presume that
$\,e^{\pm}\,$ pairs constitute the gas. Note that the pair distributions
strictly should be treated as one-dimensional due to rampant
cyclotron/synchrotron cooling perpendicular to the strong magnetic fields. Yet
addressing this nuance via the introduction of 1D Maxwellians does not
qualitatively alter the conclusions drawn here.

For relativistic pairs, the sound speed is expressed in Eq. (316) of Synge
(1957) in terms of modified Bessel functions whose argument is $\,1/\Theta\,$,
where $\,\Theta=kT/m_{e}c^{2}\,$ is the dimensionless temperature. The
observed values of $E_{\rm p}$ in the 30-35 keV range (see Figs. 1 and 2) for
the photons suggest $\,kT\sim 10\,$keV for the temperature of the photon-pair
conglomerate. This translates to $\,\Theta\sim 0.02\,$, for which the sound
speed approximately satisfies the familiar non-relativistic expectation
$\,c_{s}^{2}/c^{2}\approx 5\Theta/3\,$ (i.e. $\,\gamma P/\rho\,$ for adiabatic
index $\,\gamma=5/3\,$). It then follows that the 3D sound speed realizes
$\,c_{s}\approx 0.176c\,$. The addition of moderate radiation pressure will
increase this somewhat, but still fall well short of the relativistic EOS
result $\,c_{s}=c/\sqrt{3}\,$,

With this evaluation, the natural timescale for sound propagation over a
neutron star radius is $\,R_{\hbox{\sixrm NS}}/c_{s}\sim 0.18\,$ms for pair
plasma (with $\,R_{\hbox{\sixrm NS}}=10^{6}\,$cm). Accordingly, a $\,24\,$ms
fluctuation might signal acoustic propagation on a flux tube of length around
$\,{\cal S}\sim 130R_{\hbox{\sixrm NS}}\,$. The length might be somewhat
smaller because of the adiabatic cooling along the flux tube that lowers both
$\,\Theta\,$ and $\,c_{s}\,$ away from its footpoints at the stellar surface.
Given magnetic flux conservation along such tubes, the cross sectional area
$A$ rises at high altitudes and inversely as the field strength $B$, so that
the plasma density scales roughly as $n_{p}\propto A^{-1}\propto B$ as it
flows along a flux tube. The obvious deliverable here is that if a flux tube
with a relatively small cross sectional area constitutes the active burst
zone, then observations of $E_{\rm p}$ fluctuations can provide a measure of
the flux tube length. The arclength $\,\cal S\,$ of a dipolar field line from
footpoint to footpoint in flat spacetime geometry is stated in Eq. (23) of
Wadiasingh et al. (2018). This can easily be expressed in terms of the
footpoint colatitude $\,\theta_{\rm f}\,$, and in the domain of polar field
lines, $\,\theta_{\rm f}\ll 1\,$ yields the approximation $\,{\cal S}\approx
2.76R_{\hbox{\sixrm NS}}/\theta_{\rm f}^{2}\,$. This dipolar result would
indicate footpoint colatitudes $\,\theta_{\rm f}\sim 8^{\circ}\,$ for a
$\,24\,$ms fluctuation. This estimate increases somewhat (i.e. moving away
from the pole) when twists are introduced to modify the field morphology,
thereby lengthening field loops (profoundly near the pole: Hu et al., 2022).
This $\,\theta_{\rm f}\,$ estimate for an activated flux tube is broadly
consistent with those obtained in axisymmetric plasma simulations of twisted
magnetospheres (e.g., Chen & Beloborodov, 2017). We note that if ions are
abundant in the burst zone, then the sound speed will drop considerably,
thereby decreasing the inferred flux tube length $\,{\cal S}\,$ and moving its
locale somewhat towards the magnetic equator.

The two-blackbody emission region area estimates of $\,\sim 280-360\,$km2 for
the $\,T_{1}\sim 5\,$keV cool Planck component (see Section 3.1) can be
employed to estimate the typical cross sectional diameter $\,d\,$ of the flux
tube; the cool region most likely corresponds to the major portion of the
tube’s length $\,{\cal S}\,$ centered on its equatorial apex. Assuming that
the tube has an approximately circular cross section, its external observer-
facing surface area is of order $\,\pi{\cal S}\,d/2\,$, so that one deduces
that the typical cross sectional diameter of the tube near its apex is around
$\,1.4-1.8\,R_{\hbox{\sixrm NS}}\ll{\cal S}\,$, so that the flux tube is
appropriately slender. The corresponding tube diameter near its surface
footpoints is a small fraction of $\,R_{\hbox{\sixrm NS}}\,$ for both dipolar
and twisted field geometries.

Thus, with this intriguing and serendipitous observation, if $E_{\rm p}$
fluctuations can be confirmed in bursts, they enable a path to discerning
their emission locale. Yet the norm will be the circumstance where the active
region possesses a large number of flux tubes of different lengths and
spanning a larger cross sectional area $\,A\,$, thereby providing destructive
interference (noisiness) that obscures signature QPSO periods for acoustic
oscillations. Thus, we expect that observations of $E_{\rm p}$ oscillations
will be rare in the magnetar burst database, as appears to be the case.

##  5 Summary

In this study, we have presented the possible discovery of fluctuations of the
spectral peak energy $E_{\rm p}$ from two bursts from SGR J1935+2154 during an
outburst in early 2022. Time-integrated analysis of both bursts using
Fermi/GBM data found that their spectra were non-thermal. Time-resolved
analysis of $E_{\rm p}$ in both events provided initial evidence of $E_{\rm
p}$ oscillations (QPSOs) over the brightest part of each event, while their
lightcurve morphology presented no indication of classical QPOs. Careful
periodicity analysis supported this QPSO hypothesis through the suggestion of
a well-constrained, 24 ms (42 Hz) QPSO in one of the events (Burst 1). We
conjecture that these heretofore unprecedented events can be explained by
density and pressure perturbations propagating along a highly magnetized flux
tube that serves as an acoustic cavity. This tube is nominally of a length of
$\,{\cal S}\lesssim 130\,$ neutron star radii for a purely pair plasma in the
emission zone, though the $\,\cal S\,$ estimate drops considerably if ions are
present and thereby lower the sound speed.

Both Burst 1 and Burst 2 occurred within $\,\sim 4\,$ days of each other,
during the same outburst in January 2022. The spectroscopic and lightcurve
data shown in Fig. 1 show tantalizingly similar properties for the bursts.
This suggests, perhaps speculatively, that both events and their acoustic
oscillations arise either along the same flux tube, or a pair of flux tubes of
comparable dimensions anchored on similar magnetic field lines. This would
then provide an interesting constraint on the activation timescales over which
the environment of the magnetosphere of SGR J1935+2154  (and perhaps,
magnetars in general), change.

Detailed time-resolved analyses of other magnetar bursts are encouraged in
order to find more examples of acoustic oscillations from magnetar sources.
However, we note the difficulties in obtaining similar findings due to the
need for magnetar bursts to be sufficiently bright and long enough in order to
obtain enough oscillations or cycles, with which to confidently identify
similar or varied QPSO frequencies above the noise level. This and the
expected rarity of the proposed mechanism may make such instances rare when
searching magnetar burst databases. Additionally, we use this study to
advocate for physically accurate null hypothesis models to validate the
detection of any future QPOs and QPSOs, as the current literature and models
do not accurately describe magnetars sufficiently well.

O.J.R. gratefully acknowledges NASA funding through contract 80MSFC17M0022.
M.G.B. acknowledges the generous support of the National Aeronautics and Space
Administration through grants 80NSSC22K0777 and 80NSSC22K1576. D.H. is
supported by the Women In Science Excel (WISE) programme of the Netherlands
Organisation for Scientific Research (NWO). E.G. and Y.K. acknowledge the
support from the Scientific and Technological Research Council of Turkey
(TÜBİTAK project number 121F266).

## References

  * Anderson & Burnham (2004) Anderson, D., & Burnham, K. 2004, Second. NY: Springer-Verlag, 63, 10
  * Baring & Harding (2007) Baring, M. G., & Harding, A. K. 2007, Ap&SS, 308, 109, doi: 10.1007/s10509-007-9326-x
  * Bissaldi et al. (2009) Bissaldi, E., von Kienlin, A., Lichti, G., et al. 2009, Experimental Astronomy, 24, 47, doi: 10.1007/s10686-008-9135-4
  * Bochenek et al. (2020) Bochenek, C. D., Ravi, V., Belov, K. V., et al. 2020, Nature, 587, 59, doi: 10.1038/s41586-020-2872-x
  * Castro-Tirado et al. (2021) Castro-Tirado, A. J., Østgaard, N., Göğüş, E., et al. 2021, Nature, 600, 621, doi: 10.1038/s41586-021-04101-1
  * Chen & Beloborodov (2017) Chen, A. Y., & Beloborodov, A. M. 2017, ApJ, 844, 133, doi: 10.3847/1538-4357/aa7a57
  * CHIME/FRB Collaboration et al. (2020) CHIME/FRB Collaboration, Andersen, B. C., Bandura, K. M., et al. 2020, Nature, 587, 54, doi: 10.1038/s41586-020-2863-y
  * Hu et al. (2022) Hu, K., Baring, M. G., Harding, A. K., & Wadiasingh, Z. 2022, ApJ, 940, 91, doi: 10.3847/1538-4357/ac9611
  * Hübner et al. (2022a) Hübner, M., Huppenkothen, D., Lasky, P. D., & Inglis, A. R. 2022a, ApJS, 259, 32, doi: 10.3847/1538-4365/ac49ec
  * Hübner et al. (2022b) Hübner, M., Huppenkothen, D., Lasky, P. D., et al. 2022b, arXiv e-prints, arXiv:2205.12716.  https://arxiv.org/abs/2205.12716
  * Huppenkothen et al. (2014a) Huppenkothen, D., Heil, L. M., Watts, A. L., & Göğüş, E. 2014a, The Astrophysical Journal, 795, 114, doi: 10.1088/0004-637X/795/2/114
  * Huppenkothen et al. (2014b) Huppenkothen, D., Watts, A. L., & Levin, Y. 2014b, The Astrophysical Journal, 793, 129, doi: 10.1088/0004-637X/793/2/129
  * Huppenkothen et al. (2013) Huppenkothen, D., Watts, A. L., Uttley, P., et al. 2013, ApJ, 768, 87, doi: 10.1088/0004-637X/768/1/87
  * Huppenkothen et al. (2014c) Huppenkothen, D., D’Angelo, C., Watts, A. L., et al. 2014c, The Astrophysical Journal, 787, 128, doi: 10.1088/0004-637X/787/2/128
  * Israel et al. (2016) Israel, G. L., Esposito, P., Rea, N., et al. 2016, MNRAS, 457, 3448, doi: 10.1093/mnras/stw008
  * Kaneko et al. (2021) Kaneko, Y., Göğüş, E., Baring, M. G., et al. 2021, ApJ, 916, L7, doi: 10.3847/2041-8213/ac0fe7
  * Kass & Raftery (1995) Kass, R., & Raftery, A. 1995, Journal of the American Statistical Association, 90, 773
  * Kirsten et al. (2021) Kirsten, F., Snelders, M. P., Jenkins, M., et al. 2021, Nature Astronomy, 5, 414, doi: 10.1038/s41550-020-01246-3
  * Kozyrev et al. (2022) Kozyrev, A. S., Golovin, D. V., Litvak, M. L., et al. 2022, GRB Coordinates Network, 31512, 1
  * Levin (2007) Levin, Y. 2007, MNRAS, 377, 159, doi: 10.1111/j.1365-2966.2007.11582.x
  * Li et al. (2022) Li, X., Ge, M., Lin, L., et al. 2022, The Astrophysical Journal, 931, 56, doi: 10.3847/1538-4357/ac6587
  * Liddle (2007) Liddle, A. R. 2007, Monthly Notices of the Royal Astronomical Society: Letters, 377, L74, doi: 10.1111/j.1745-3933.2007.00306.x
  * Lin et al. (2020) Lin, L., Göğüs, E., Roberts, O. J., et al. 2020, ApJ, 902, L43, doi: 10.3847/2041-8213/abbefe
  * Lin et al. (2011) Lin, L., Kouveliotou, C., Baring, M. G., et al. 2011, The Astrophysical Journal, 739, 87, doi: 10.1088/0004-637X/739/2/87
  * Lin et al. (2012) Lin, L., Göğüş, E., Baring, M. G., et al. 2012, The Astrophysical Journal, 756, 54, doi: 10.1088/0004-637X/756/1/54
  * Meegan et al. (2009) Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791, doi: 10.1088/0004-637X/702/1/791
  * Mereghetti et al. (2020) Mereghetti, S., Savchenko, V., Ferrigno, C., et al. 2020, ApJ, 898, L29, doi: 10.3847/2041-8213/aba2cf
  * Miller et al. (2019) Miller, M. C., Chirenti, C., & Strohmayer, T. E. 2019, The Astrophysical Journal, 871, 95, doi: 10.3847/1538-4357/aaf5ce
  * Pavlović et al. (2013) Pavlović, M. Z., Urošević, D., Vukotić, B., Arbutina, B., & Göker, Ü. D. 2013, ApJS, 204, 4, doi: 10.1088/0067-0049/204/1/4
  * Ridnaia et al. (2022) Ridnaia, A., Frederiks, D., Ulanov, M., et al. 2022, GRB Coordinates Network, 31497, 1
  * Ridnaia et al. (2021) Ridnaia, A., Svinkin, D., Frederiks, D., et al. 2021, Nature Astronomy, 5, 372, doi: 10.1038/s41550-020-01265-0
  * Roberts (2023) Roberts, O. J. 2023, In prep.
  * Roberts & Fermi GBM Team (2022) Roberts, O. J., & Fermi GBM Team. 2022, GRB Coordinates Network, 31445, 1
  * Roberts et al. (2021) Roberts, O. J., Veres, P., Baring, M. G., et al. 2021, Nature, 589, 207, doi: 10.1038/s41586-020-03077-8
  * Schwarz (1978) Schwarz, G. 1978, The Annals of Statistics, 6, 461
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stamatikos et al. (2014) Stamatikos, M., Malesani, D., Page, K. L., & Sakamoto, T. 2014, GRB Coordinates Network, 16520, 1
  * Strohmayer & Watts (2005) Strohmayer, T. E., & Watts, A. L. 2005, The Astrophysical Journal, 632, L111, doi: 10.1086/497911
  * Strohmayer & Watts (2006) —. 2006, The Astrophysical Journal, 653, 593, doi: 10.1086/508703
  * Synge (1957) Synge, J. L. 1957, The Relativistic Gas
  * von Kienlin et al. (2020) von Kienlin, A., Meegan, C. A., Paciesas, W. S., et al. 2020, ApJ, 893, 46, doi: 10.3847/1538-4357/ab7a18
  * Wadiasingh et al. (2018) Wadiasingh, Z., Baring, M. G., Gonthier, P. L., & Harding, A. K. 2018, ApJ, 854, 98, doi: 10.3847/1538-4357/aaa460
  * Younes et al. (2017) Younes, G., Kouveliotou, C., Jaodand, A., et al. 2017, The Astrophysical Journal, 847, 85, doi: 10.3847/1538-4357/aa899a
  * Younes et al. (2020) Younes, G., Güver, T., Kouveliotou, C., et al. 2020, The Astrophysical Journal Letters, 904, L21, doi: 10.3847/2041-8213/abc94c
  * Younes et al. (2022) Younes, G., Baring, M. G., Harding, A. K., et al. 2022, arXiv e-prints, arXiv:2210.11518.  https://arxiv.org/abs/2210.11518
  * Zhong et al. (2020) Zhong, S.-Q., Dai, Z.-G., Zhang, H.-M., & Deng, C.-M. 2020, ApJ, 898, L5, doi: 10.3847/2041-8213/aba262


# Some properties of Euler capital allocation

Lars Holden  Norwegian Computing Center, P. O. Box 114 Blindern, NO–0314 Oslo,
Norway<EMAIL_ADDRESS>
(Date: )

###### Abstract.

The paper discusses capital allocation using the Euler formula and focuses on
the risk measures Value-at-Risk (VaR) and Expected shortfall (ES). Some new
results connected to this capital allocation is shown. Two examples illustrate
that capital allocation with VaR is not monotonous which may be surprising
since VaR is monotonous. A third example illustrates why the same risk measure
should be used in capital allocation as in the evaluation of the total
portfolio. We show how simulation may be used in order to estimate the
expected Return on risk adjusted capital in the commitment period of an asset.
Finally, we show how Markov chain Monte Carlo may be used in the estimation of
the capital allocation.

##  1\. Introduction

The regulatory framework requires a quantification of the total risk in a
corporate. Based on the quantification and a risk measure there is a
requirement for an economic capital that is able to absorb potential losses.
The competition makes it necessary for a financial corporate to ensure an
efficient use of their economic capital. As a part of improving the use of the
economic capital, the capital that is allocated to each asset is calculated.
This makes it possible to evaluate each part of the portfolio. This paper
studies some properties of capital allocation assuming that the stochastic
properties of the different assets including their correlations are known.
Hence, we will not discuss how to estimate these distributions.

The two papers Artzner et al. [2] and [3] have initiated a large number of
papers giving a much better understanding of capital allocation. There seem to
be an agreement that capital should be allocated proportional with the partial
derivatives of the risk measure and use the Euler Theorem, see Tasche [12].
This was first proposed by Litterman [7]. Denault [4] argues for the same
capital allocation using theory from cooperative games. Kalkbrener et al., [6]
shows how axioms assuming that the capital allocation is linear and
diversifying lead to the same capital allocation. Tasche [13] and [14] and
Fischer [5] give an overview of the argument for using Euler allocation.

A large number of authors, see e.g. Acerbie et al [1], argue that the risk
measure should be coherent. This includes for example expected shortfall (ES)
but excludes Value-at-Risk (VaR), the most used risk measure. VaR is used as a
part of the Basel II framework. VaR is not subadditive implying that there may
be an additional cost when adding two portfolios instead of a saving due to
diversification. Several authors argue that one should use ES which is the
smallest measure that is law-permitting dominating VaR, see Tasche [12].

We will contribute in this area with several different examples. First
illustrating that capital allocation using VaR is not monotonous. This may be
surprising since VaR is monotonous. It may be a larger problem than the
missing subadditivity of VaR in a regular use of VaR in capital allocation in
corporate. This easily leads to suboptimal performance of the management of
the portfolio.

The regulatory framework may require that VaR is used for estimating the
economic capital needed for a portfolio. It may however be tempting to use ES
for the capital allocation of each asset since this is a coherent measure with
better mathematical properties. We give an example illustrating that the same
risk measure should be used for the total portfolio and the capital
allocation, else one easily gets conflict of interest.

Each part of the portfolio may be a commitment for a longer period, and the
corporate will not always be able to end the part of the portfolio that is not
cost efficient based on the latest evaluations. Hence, a new investment should
be evaluated based on the expected capital requirement in the commitment
period for the investment, not only based on today’s portfolio. Capital
allocation based on the present portfolio is not optimal since the total
portfolio may change during the period where it is not possible to change the
condition for a part of the portfolio. The natural alternative is to simulate
the portfolio and use the expected average value of the economic capital for
the commitment period of the new investment. This is illustrated in an
example.

The numerical calculations involved in a capital allocation may be very
computer intensive. Direct Monte Carlo simulation may not work since we need a
large number realizations to make a good evaluation of the tail behavior.
Importance sampling reduces the problem, see Kalkbrener et al. [6], but it may
still be too CPU requiring. We propose to use Markov chain Monte Carlo. Then
we may focus such that all realizations are used in the quantification. This
is illustrated in an example.

##  2\. Model

In the following we use the terminology from Tasche [14]. Let the random
variable $X_{i}$ denote the cash flow from asset $i$. A portfolio $X$ consists
of $n$ different assets

$X=\sum_{i=1}^{n}X_{i}.$

Further, let $EC(X)=\rho(X)$ denote the economic capital that is deemed
necessary by the regulator or corporate in order to handle a possible
unfortunate development of the portfolio. $\rho(X)$ is a risk measure on the
portfolio $X\in V$ where $V$ is the set of real valued random variables. We
may define $\rho$ as VaR, ES, the standard deviation multiplied with a
constant or any other measure for the uncertainty. We need to define some
properties on risk measures. A risk measure is monotonous if $X,Y\in V$ and
$X\leq Y$ almost everywhere implies

$\rho(X)\geq\rho(Y),$

subadditive  if $X,Y\in V$ implies

$\rho(X+Y)\leq\rho(X)+\rho(Y),$

positive homogeneous if $X\in V$, $h\in R$, and $h>0$ implies

$\rho(hX)=h\rho(X)$

and translation invariant if $X\in V$ and $h\in R$ implies

$\rho(X+h)=\rho(X)-h.$

A risk measure is denoted coherent  if it is monotonous, subadditive, positive
homogeneous and translation invariant. In this paper we will assume the risk
measure is positive homogeneous but do not require the other properties above.

In order to define the derivative of $\rho$ we introduce the function

$f({\bf u})=\rho(\sum_{i=1}^{n}u_{i}X_{i})$

where ${\bf u}=(u_{1},\cdots,u_{n})$. Assuming the risk measure is positive
homogeneous, it satisfies the Euler formula, see Tasche [12]

$f({\bf u})=\sum_{i=1}^{d}u_{i}\frac{\partial f}{\partial u_{i}}({\bf u})$

where we assume the partial derivatives exists. This makes it natural to
allocate the capital

$\rho(X_{i},X)=\frac{\partial f}{\partial u_{i}}({\bf u}_{1})$

where ${\bf u}_{1}=(1,\cdots,1)$ to the asset $X_{i}$. This capital allocation
is denoted Euler allocation. We see immediately that this capital allocation
is linear, e.g.

$\rho(X_{i}+X_{j},X)=\rho(X_{i},X)+\rho(X_{j},X)$

and has the full allocation property, e.g.

$\rho(X)=\sum_{i=1}^{n}\rho(X_{i},X).$

Assuming $\rho(X)$ is coherent, Denault [4] has proved using game theory, that
Euler allocation is the only allocation satisfying that we always have
$\rho(X_{i},X)\leq\rho(X_{i})$. Kalkbrener et al., [6], has proved a similar
result using axioms. This is a critical property closely connected to
subadditivity. It states that an asset is allocated less capital as part of a
portfolio than alone. Tasche [11] and [14] give an overview of the argument
for this capital allocation.

Another measure on the portfolio is the Return on risk adjusted capital,

$\mbox{RORAC}(X)=\frac{EX}{EC(X)}.$

For each asset $X_{i}$ in a portfolio $X$ we define
$\mbox{RORAC}(X_{i},X)=EX_{i}/\rho(X_{i},X).$ A capital allocation is denoted
RORAC compatible, see Tasche [11], if there exists $\varepsilon_{i}>0$ such
that $\mbox{RORAC}(X_{i},X)>\mbox{RORAC}(X)$ implies

$\mbox{RORAC}(X+hX_{i})>\mbox{RORAC}(X)$

for all $0<h<\varepsilon_{i}$. The intuition is that if asset $X_{i}$ has
higher RORAC than the entire portfolio, then increasing this asset in the
portfolio increases the RORAC of the portfolio. This is illustrated in Example
3.

Tasche [9] proves that Euler allocation is the only allocation that is RORAC
compatible. In the following two sections we will discuss Euler allocation for
two of the most popular risk measures.

##  3\. Value-at-Risk, VaR

Defined the $\alpha$-quantile $q_{\alpha}$ as

$q_{\alpha}(X)=\inf\\{z\in R|P(X\leq z)\geq\alpha\\}.$

Then the risk measure VaR is defined as $VaR_{\alpha}(X)=q_{\alpha}(-X)$. We
typically have $\alpha\geq 0.99$ implying that estimation of VaR only focus on
one point in the tail. This makes estimates for VaR much less stable than for
example standard deviation.

VaR is monotone, positive homogeneous, and translation invariant. But it is
also well-known that VaR is not subadditive since we may have

$VaR_{\alpha}(X+Y)>VaR_{\alpha}(X)+VaR_{\alpha}(Y).$

This implies that VaR is not a coherent risk measure. VaR expresses the
economical capital necessary to ensure that the probability for a default is
less than $\alpha$. Hence, VaR focuses only on one point in the cumulative
distribution of $X$, the maximum value $z$ where $P(X\geq z)\geq\alpha$. This
property is not additive. This implies that if we have two portfolios $X$ and
$Y$ the capital requirement may be larger when we add them together than if we
keep them separate. This may be an argument for splitting the portfolio and
hence the corporate in two. This may seem counter-intuitive as a risk measure
and a not wanted property for the corporate for strategic reasons. There are
several papers stating that risk measures that are not coherent should not be
used. Artzner et al. [3] give three examples with discrete random variables,
one illustrating that the measure is not subadditive, another illustrating
that VaR fails to recognize concentration of risk and fails to encourage a
reasonable allocation of risk between agents. Tasche [10] gives an example
with two independent Pareto distributions that does not satisfy subadditivity.
Acerbi et al, [1], writes “ ..if a measure is not coherent we just choose not
to call it a risk measure at all”, particular due to the missing subadditive
property and since there exists coherent risk measures with satisfactory
properties.

The Euler allocation for VaR is

(1) $VaR_{\alpha}(X_{i},X)=-E\\{X_{i}|X=-VaR_{\alpha}(X)\\}$

under some smoothness assumptions, see Tasche [9]. Kalkbrener et al. [6]
report that capital allocation with VaR may require larger economic capital to
an asset then the lowest possible outcome of the variable. But then the
capital allocation is using covariance instead of Euler allocation. This will
not happen with Euler capital allocation as is seen from (1). It is much
easier to perform capital allocation using correlation than Euler allocation.
In Section 7 we show how calculation of capital allocation using Euler
allocation may be performed efficiently with Markov chain Monte Carlo methods.

VaR satisfies the following monotonicity property: If $P(X\leq z)\geq P(Y\leq
z)$ for all values of $z\in R$ then $VaR_{\alpha}(X)\geq VaR_{\alpha}(Y)$.
However, Euler capital allocation with VaR does not satisfy the same
monotonicity e.g. we may have $VaR_{\alpha}(X,X+Y)\leq VaR_{\alpha}(Y,X+Y)$.
This may be surprising. Missing this property may be more important than
missing subadditivity when considering whether the capital allocation inside a
corporate is fair or not. To the authors knowledge, this property is not
proved earlier. This property is due to the fact that VaR focuses on only one
point in the distribution.

We give two examples that Euler capital allocation with VaR is not monotonous.
The first example has two independent discrete variables and the second
example has stochastic variables that are continuous and dependent.

Example 1  Let $X_{1}$ and $X_{2}$ be two independent stochastic variables
where

$P(X_{1}=0)=P(X_{2}=0)=0.9925$

and

$P(X_{1}=-200)=P(X_{2}=-100)=0.0075.$

Then

$\displaystyle VaR_{0.99}(X_{1})$ $\displaystyle=VaR_{0.99}(X_{2})=0$
$\displaystyle VaR_{0.99}(X_{1}+X_{2})$ $\displaystyle=100$ $\displaystyle
VaR_{0.99}(X_{1},X_{1}+X_{2})$ $\displaystyle=0$ $\displaystyle
VaR_{0.99}(X_{2},X_{1}+X_{2})$ $\displaystyle=100$

Hence, the asset $X_{2}$ gets allocated all the risk even though $X_{1}\leq
X_{2}$. This may seem surprising and the example is analysed more thoroughly
by introducing weights $(u_{1},u_{2})$ for the two assets, e.g.

$X(u_{1},u_{2})=u_{1}X_{1}+u_{2}X_{2}.$

Then we have

$\displaystyle VaR_{0.99}(X(u_{1},1))$ $\displaystyle=100u_{1}\mbox{ \ \ for
$0\geq u_{1}\geq 2$ }$ $\displaystyle VaR_{0.99}(X(1,u_{2}))$
$\displaystyle=100\mbox{ \ \ for $u_{2}\geq 0.5.$ }$

We see that the risk for the portfolio is sensitive to weight $u_{1}$ but not
$u_{2}$ in the interesting region close to $(u_{1},u_{2})=(1,1)$, hence it
makes sense to reduce the weight of $X_{1}$ instead of $X_{2}$ if the
regulatory requirement is based on $VaR_{0.99}(X).$ This also illustrates that
it is possible to cheat in the system. The person responsible for asset
$X_{1}$ may reduce the capital allocated to this asset if he promises to give
a value of 100 to charity if $X_{1}=-100.$ This makes the distribution of
$X_{1}$ and $X_{2}$ equal implying that they get the same capital allocation.

Example 2  Let $X=3X_{1}+X_{2}$ where $X_{1}$ is symmetric around $0$ and
$X_{2}$ is given as

$X_{2}=\left\\{\begin{array}[]{cc}-X_{1}&\mbox{ if $X_{1}\leq 0$ }\\\
-2X_{1}&\mbox{ if $X_{1}>0$. }\end{array}\right.$

This implies that $X_{2}$ has the same upside as $X_{1}$, but the downside is
not as good considered as a univariate variable. But $X_{2}$ has good
diversification properties in the portfolio. The portfolio $X$ has the same
marginal distribution as $X_{2}$,

$X=\left\\{\begin{array}[]{cc}2X_{1}&\mbox{ if $X_{1}\leq 0$ }\\\ X_{1}&\mbox{
if $X_{1}>0$. }\end{array}\right.$

Then $VaR_{\alpha}(X_{2},X)<0<VaR_{\alpha}(X_{1},X)$ even though $P(X_{1}\leq
z)\leq P(X_{2}\leq z)$ all $z\in R$.

The situation is very different in the two cases. In Example 1 it is tempting
to deviate from Euler capital allocation by moving capital between the
involved assets in order to maintain monotonicity. In Example 2 it is
important to not change the Euler capital allocation in order to honor the
increased diversification.

##  4\. Expected shortfall, ES

ES is defined as

$ES_{\alpha}(X)=-E\\{X|X\leq-VaR_{\alpha}(X)\\}.$

Acerbi et al, [1] prove that

$ES_{\alpha}(X)=\frac{1}{1-\alpha}\int_{\alpha}^{1}VaR_{\tau}(X)d\tau.$

ES is a coherent risk measure. Further, it is proved that Euler allocation
gives

$ES_{\alpha}(X_{i},X)=-E\\{X_{i}|X\leq-VaR_{\alpha}(X)\\}$

assuming sufficient smoothness. Since ES takes the average of VaR for
$\alpha\leq\tau\leq 1$ we have $ES_{\alpha}(X)\leq VaR_{\alpha}(X)$. It is
also proved the ES is the smallest risk measure that is law invariant where
$\rho(X)\leq VaR_{\alpha}(X)$. This makes ES a good alternative to VaR. The
property law invariant  is that if $X,Y\in V$ and $P(X\leq z)=P(Y\leq z)$ for
all $z\in R$ then $\rho(X)=\rho(Y).$ Since ES focuses on the entire tail and
not only a quantile, it may require more realizations in a Monte Carlo
estimation than VaR.

##  5\. Combining ES and VaR

Since ES is coherent while VaR is not, many authors recommend to use ES
instead of VaR. We will discuss under the assumption that the regulator
requires the use of VaR for setting the economic capital whether to use VaR or
ES in the capital allocation. Since VaR is not subadditive some authors seem
to recommend to use ES for capital allocation. The missing monotonicity
property strengthen this argument. However, we will argue that this is not
necessarily a good choice. It is correct that VaR may lead to examples where
we do not get subadditivity and therefore there may be arguments for splitting
the portfolio in several parts. However, most practitioners report a 30%
diversification effect indicating that the missing subadditivity property is
mainly academic. The following example shows that combining VaR as a risk
measure for the portfolio and ES in the capital allocation, may give conflict
of interest. This is avoided if VaR or ES is used both as a risk measure for
the portfolio and in the capital allocation. This indicates that if VaR is
used as risk measure, then VaR should also be used in the capital allocation.

Figure 1.  RORAC for the porfolio (line) and the assets $uX_{1}$ (dotted) and
$(1-u)X_{1}$ (dashed) with Euler capital allocation using VaR.  Figure 2.
RORAC for the porfolio (line) and the assets $uX_{1}$ (dotted) and
$(1-u)X_{1}$ (dashed) with Euler capital allocation using ES. Figure 3.  RORAC
for the porfolio (line) and the assets $uX_{1}$ (dotted) and $(1-u)X_{1}$
(dashed) with VaR as risk measure and $\rho_{VaR-ES}(X_{i},X)$ used in the
capital allocation.

Example 3  Let $X_{i}=0.5-I_{i}Y_{i}$ for $i=1,2$ where $I_{i}$ is an
indicator and $Y_{i}$ is a Pareto distributed variable with density

$f(y)=\frac{\gamma_{i}}{b_{i}}(\frac{y}{b_{i}}+1)^{-\gamma_{i}-1}$

where $y>0.$ The indicator $I_{i}=1$ with probability $0.1$ and else
$I_{i}=0$. We have chosen the variables $\gamma_{1}=5$ and $\gamma_{2}=1.7$
and $b_{i}$ such that $E\\{X_{i}\\}=0.2$ for $i=1,2$. We study the portfolio
$X=uX_{1}+(1-u)X_{2}$ for $0\leq u\leq 1$ and compare RORAC with Euler capital
allocation using VaR and ES as risk measure (Figures 1 and 2), and when VaR is
used as risk measure for the portfolio while

$\rho_{VaR-ES}(X_{i},X)=ES_{\alpha}(X_{i},X)VaR_{\alpha}(X)/ES_{\alpha}(X)$

is used in the capital allocation ( Figure 3). The figures show RORAC for the
portfolio and the two assets $uX_{1}$ and $(1-u)X_{1}$ as a function of $u$
for the different risk measures and capital allocation methods. When we use
Euler capital allocation with VaR or ES we get an optimal portfolio for
$u=0.3$ and $u=0.7$ respectively where both assets and the total portfolio
have the same RORAC. That the assets have the same RORAC in the optimum as the
total portfolio follows from the RORAC compatibility of Euler capital
allocation. In these two cases RORAC is increasing for the assets when the
asset gets less weight in the portfolio. When $\rho_{VaR-ES}$ is used for the
capital allocation the performance is more complex. The diversification effect
is not as expected and the two assets have very different RORAC at the optimal
diversification. If we require that the two assets have approximately the same
RORAC then the total portfolio is far from optimal. If we reduce the weight of
the asset with lowest RORAC, we may reduces the RORAC of the total portfolio.

##  6\. Changes in the portfolio

Assume we already have a portfolio $X$ and consider to expand the portfolio
with $X_{n+1}$. Then we should base the evaluation on the capital allocation

(2) $\rho(X_{n+1},X+X_{n+1}).$

Euler capital allocation does not encourage new investments $X_{n+1}$ where

$\rho(X,X+X_{n+1})-\rho(X)$

is large, hence contributes to a reduction in the capital allocated to the
rest of the portfolio. There are capital allocation methods that takes these
other properties into account. But Tasche [12] proves that for continuous
differentiable, subadditive and positive homogeneous risk measures then Euler
allocation gives

$\rho(X_{n+1},X+X_{n+1})\geq\rho(X+X_{n+1})-\rho(X).$

Hence if each asset only gets allocated the additional increase in the total
portfolio, then the sum of the capital allocated to each asset may not add up
to the capital required for the entire portfolio and the sum will never be
above. In some cases it may be fair that the last asset, $X_{n+1}$, gets
allocated only the marginal increase and all the other assets get allocated
capital assuming that $X_{n+1}$ is not present. But then it is necessary to
argue that this particular asset is different from the rest and deserves a
special treatment.

In a one-period framework, as pointed out by Fischer [5], where we at time
$t=0$ can evaluate the entire portfolio $X$ and we are not allowed to do any
changes before $t=T$ and the variables $X$ and $X_{n+1}$ denote the values at
time $t=T$, then it is natural to base the evaluation on (2). Euler capital
allocation is only based on the risk as it is evaluated today. The typical
situation is often very different. It does not consider that a large part of
the portfolio may be invested a long time ago, under different economic
conditions and that the corporate may not be in the position to remove part of
the portfolio that is not beneficial any more. The present decision is whether
to include a new investment in the portfolio. Let $X_{t}$ and $X_{n+1,t}$
denote the portfolio and the new investment at time $t$ and assume the new
investment is a commitment from $t=0$ to $t=T.$ Then it is natural to base the
decision on the capital allocation

$\frac{1}{T}\int_{t=0}^{T}E\\{\rho(X_{n+1,t},X_{t}+X_{n+1,t})\\}dt.$

In order to evaluate this expression it is necessary to evaluate the future
properties of the portfolio and the new investment. This is typically done by
simulation. The following schematic example shows how this may be done.

Table 1.  RORAC for a new investment depending on the other investment. New investment | Other investment
---|---
| $X_{1}$ | $X_{2}$
$X_{1}$ | 0.030 | 0.033
$X_{2}$ | 0.071 | 0.044
Table 2.  Average RORAC for a two years period for a new investment depending on the old investment in the portfolio New investment | Old investment
---|---
| $X_{1}$ | $X_{2}$
$X_{1}$ | 0.031 | 0.032
$X_{2}$ | 0.064 | 0.051

Example 4  Let the portfolio be $X=Z_{1}+Z_{2}$ where one of the $Z_{j}$ is
renewed each year for a period of two years. Each of the variables $Z_{j}$ may
be of the two types $X_{i}$, $i=1,2$ as defined in Example 3\. When we make
the investment for a two years period we know the portfolio for the first
year, but for the second year the other part of the portfolio is with equal
probability equal to $X_{1}$ or $X_{2}$. Table 1 shows RORAC for the four
cases. If the new investment is of type $X_{2}$ and the other investment that
is lasting for one more year is also of type $X_{2}$, then the RORAC for the
new investment the first year is 0.044. For the second year will the other
investment be of type $X_{1}$ or $X_{2}$ with equal probability leading to a
RORAC equal 0.071 or 0.044 with equal probability. This implies that the
expected average annual RORAC for the new investment is equal to
(0.044+(0.044+0.071)/2)/2=0.51. By a similar method we may calculate the
numbers shown in Table 2. This table may be used in the evaluation of a new
investment based on whether the new and the old investment is of type $X_{1}$
or $X_{2}$.

##  7\. Calculation of capital allocation

The easiest method to calculate the Euler capital allocation based on the
stochastic properties of the portfolio and a risk measure is to use Monte
Carlo simulation. Then it is generated a large number of realizations and the
risk measure is calculated by evaluating the realizations empirically. If VaR
or ES is used as risk measure with a high value $\alpha$, we will only use a
small part of the simulated values. This makes the simulation inefficient.
This may be considerably improved by using importance sampling, see e.g.
Kalkbrener et al. [6] or Tasche [13]. If we use Markov chain Monte Carlo
(MCMC) it is possible to improve the sampling even more. See Meyn and Tweedie
[8] for a thorough introduction to MCMC. In the following example we
illustrate how capital allocation may be calculated using Monte Carlo
simulation, Importance sampling and MCMC. We give complete algorithms for each
type. It is outside the scope of this paper to discuss the different
alternatives within each class in detail.

Example 5  Let

$X^{j}=\sum_{i=1}^{n}X^{j}_{i}$

denote realization number $j$ of the portfolio. Assume

$X^{j}_{i}=a_{i}-exp(Y^{j}_{i})$

where $Y^{j}_{i}\sim N(\mu_{i},\sigma_{i}^{2})$ and all $Y^{j}_{i}$ are
uncorrelated in order to simplify the notation. Let $X^{(j)}$ denote the
$X^{j}$ sorted such that

$X^{(1)}\leq X^{(2)}\leq\cdots\leq X^{(n)}$

and $X^{(j)}_{i}$ is the realization for asset $i$ that correspond to
$X^{(j)}$. Assume $n_{\alpha}=(1-\alpha)n$ is an integer such that
$-X^{(n_{\alpha})}$ is a reasonable estimator for $VaR_{\alpha}(X)$. Let $b$
be an integer such that $b\leq n_{\alpha}$.

The Monte Carlo simulation algorithm below gives n realizations of $X$. From
these realizations we find the estimators $VaR_{\alpha,MC}(X)$ for
$VaR_{\alpha}(X)$ and $VaR_{\alpha,MC}(X_{i},X)$ for $VaR_{\alpha}(X_{i},X)$.

  1. (1)

For $j=1,2,\cdots,m$

    1. (a)

Set $X^{j}=\sum_{i=1}^{n}(a_{i}-exp(Y^{j}_{i}))$

  2. (2)

Sort to find $X^{(1)},\cdots,X^{(m)},$ and the corresponding $X^{(j)}_{i}$

  3. (3)

Set $VaR_{\alpha,MC}(X)=-X^{(n_{\alpha})}$

  4. (4)

Set
$VaR_{\alpha,MC}(X_{i},X)=-\frac{1}{2b+1}X^{(n_{\alpha})}\sum_{j=-b}^{b}\frac{X^{(n_{\alpha}+j)}_{i}}{X^{(n_{\alpha}+j})}$

In order to make good estimates for the capital allocation, it is necessary
with many realizations estimating this, hence $b$ large. On the other hand it
is necessary with $n$ large such that the $2b+1$ realizations
$X^{(n_{\alpha}-b)},\cdots,X^{(n_{\alpha}+b)}$ all are sufficiently close to
the quantile $VaR_{\alpha}(X)$. The storage problem may be solved by only
storing realizations in an interval $d_{1}\leq X^{j}\leq d_{2}$ surrounding
$VaR_{\alpha}(X)$.

Only a small portion of the realizations is in the interesting interval. This
may be improved by importance sampling. Here is a very simply importance
sampling algorithm where $Y^{j}_{i,IS}\sim N(\mu_{i,IS},\sigma_{i,IS}^{2})$.
Define $\phi_{\mu,\sigma^{2}}(y)$ as the density in the $N(\mu,\sigma^{2})$
distribution and $p^{j}$ the importance sampling weight for each realization.
Further, define $p^{(j)}$ as the values of $p^{j}$ sorted according to the
size of $X^{j}$ and scaled such that $0=p^{(0)}<p^{(1)}<\cdots<p^{(m)}=1$ and
$p^{(j+1)}-p^{(j)}=cp^{k}$ where $p^{k}$ is the corresponding value before the
sorting. Define $n_{\alpha,IS}$ such that $p^{(n_{\alpha,IS})}$ is closest to
$(1-\alpha)n$. We may then define the importance sampler estimators
$VaR_{\alpha,IS}(X)$ and $VaR_{\alpha,IS}(X_{i},X)$ by the algorithm:

  1. (1)

For $j=1,2,\cdots,m$

    1. (a)

Set $X^{j}=\sum_{i=1}^{n}(a_{i}-exp(Y^{j}_{i}))$

    2. (b)

Set
$p^{j}=\prod_{i=1}^{n}\frac{\phi_{\mu_{i},\sigma_{i}^{2}}(Y^{i}_{j})}{\phi_{\mu_{i,IS},\sigma_{i,IS}^{2}}(Y^{i}_{j})}$

  2. (2)

Sort to find $X^{(1)},\cdots,X^{(m)},$ and the corresponding $X^{(j)}_{i}$ and
$p^{(j)}$.

  3. (3)

Set $VaR_{\alpha,MC}(X)=-X^{(n_{\alpha,IS})}$.

  4. (4)

Set
$VaR_{\alpha,MC}(X_{i},X)=-\frac{1}{2b+1}X^{(n_{\alpha,IS})}\sum_{j=-b_{IS}}^{b_{IS}}\frac{X^{(n_{\alpha,IS}+j)}_{i}}{X^{(n_{\alpha}+j})}$.

Still, only a small part of the realizations are used in the estimation and
none of these are exactly such that $X=VaR_{\alpha}(X).$ By using Markov chain
Monte Carlo MCMC it is possible to make all realizations equal to
$X=VaR_{\alpha}(X)$ and such that all realizations may be used in the capital
allocation. Assume $VaR_{\alpha}(X)$ is found by another algorithm, for
example one of the two algorithms described above. It is necessary with
significant larger $n$ in order to find the capital allocation than VaR hence
it make sense to use on of the other algorithms for this. Then the MCMC
estimator for the capital allocation $VaR_{\alpha,MCMC}(X_{i},X)$ may be found
by the following algorithm:

  1. (1)

Let $X^{0}$ be any realization satisfying $X^{0}=VaR_{\alpha}(X)$.

  2. (2)

Let $k_{1}\leq n$ be any index

  3. (3)

For $j=1,2,\cdots,m*n$

    1. (a)

Find index $k_{2}\leq n$ and $k_{1}\neq k_{2}$ uniformly

    2. (b)

Find $\tilde{Y}^{j+1}_{k_{2}}\sim N(\mu_{k_{2}},\sigma_{k_{2}}^{2})$ and
$corr(\tilde{Y}^{j+1}_{k_{2}},Y^{j}_{k_{2}})$ constant

    3. (c)

Set $\tilde{X}^{j}_{k_{2}}=(a_{k_{2}}-exp(Y^{j}_{k_{2}}))$

    4. (d)

Set $\tilde{X}^{j}_{k_{1}}=X^{j}_{k_{1}}+X^{j}_{k_{2}}-\tilde{X}^{j}_{k_{2}}$

    5. (e)

Set $\tilde{Y}^{j}_{k_{1}}=\log(a_{k_{1}}-\tilde{X}^{j}_{k_{1}})$

    6. (f)

Set
$p_{j}=\phi_{\mu_{i},\sigma_{i}^{2}}(\tilde{Y}^{i}_{j})/\phi_{\mu_{i},\sigma_{i}^{2}}(Y^{i}_{j})$

    7. (g)

Set $X^{j+1}=\tilde{X}^{j+1}$ with probability $\min\\{1,p_{j}\\}$ and else
set $X^{j+1}=X^{j}$.

    8. (h)

Set $k_{1}=k_{2}$

  4. (4)

Set $VaR_{\alpha,MCMC}(X_{i},X)=-\frac{1}{m}\sum_{j=1}^{m}X^{jn}_{i}$.

The MCMC algorithm generates a chain of realizations $X^{j}$ that are
according to the distribution we want to study. But the first realizations in
the chain are not from the distribution (a burn-in period) and realizations
$X^{j}$ and $X^{j+k}$ are dependent but the dependence decreases geometrically
in $k$. Table 3 shows the number of numerical operations for the different
algorithms. Notice that the number of operations is considerably smaller for
MCMC in the capital allocation per realization.

Table 3.  Number of numerical operations (additions, multiplications, exponential, logarithm ) per realization in the three algorithms  | MC | IS | MCMC
---|---|---|---
$VaR(X)$ | 3n | 6n | -
$VaR(X_{i},X)$ | $3n/(2b+1)$ | $6n/(2b_{IS}+1)$ | 9

We have tested the three algorithm in an example with $n=90$ assets with the
following parameters: $\sigma_{i}=0.5$ and $E\\{X_{i}\\}=0.2$ for all $i$
while the assets are divided into three groups with 30 assets in each group
where $\mu_{i}=0.44/0.45/0.47$ respectively. The parameters $a_{i}$ are
determined from the other parameters. Table 4 shows the estimate for
$VaR_{0.99}(X)$ and the capital allocation $VaR_{0.99}(X_{i},X)$ for the
portfolio for the three different algorithms. By having 30 equal assets in the
portfolio, it is possible to estimate the precision of the different
estimates.

Table 4.  Estimate for the VaR and the capital allocation using the three algorithms. The standard deviation of the estimates are given en parenthesis.  | MC | IS | MCMC
---|---|---|---
m | 1.000.000 | 1.000.000 | 100.000
b | 1.600 | 20.000 | -
$VaR_{0.99}(X)$ | 6.33 | 6.38 | -
$VaR_{0.99}(X_{1},X)$ | 0.038 (0.018) | 0.042 (0.016) | 0.038 (0.027)
$VaR_{0.99}(X_{31},X)$ | 0.064 (0.021) | 0.065 (0.016) | 0.066 (0.029)
$VaR_{0.99}(X_{61},X)$ | 0.109 (0.019) | 0.108 (0.021) | 0.109 (0.026)

In importance sampling we have used $\mu_{i,IS}=\mu_{i}+0.2$ resulting in
about 10 times as many realizations in an interval close to $VaR_{0.99}(X)$ as
in Monte Carlo simulation. This implies that importance sampling may give an
improvement with a factor 10 compared to Monte Carlo simulation. In MCMC we
have used correlation equal to 0.3 when proposing a new
$\tilde{Y}^{j+1}_{k_{2}}$ value leading to an acceptance ratio at 0.57 and
that realizations with $X^{j}$ and $X^{j+5}$ are almost independent. There is
almost no burn-in in this MCMC algorithm. Importance sampling gives 2b=40.000
samples from m= 1 mill. MCMC gives 100.000 samples using about 10 % of the
number of operations. This indicates that that MCMC may give an improvement
with a factor 100 compared to Monte Carlo simulation. In the test, we have not
performed any optimization of the algorithm and parameters. Hence, it is
reasonable to assume that it is possible to improve both importance sampling
and MCMC even further.

##  8\. Concluding remarks

In this paper we have given an overview over capital allocation using the
Euler formula. We have contributed to the research showing that:

  * •

Euler capital allocation with VaR is not monotone even though VaR is monotone.
This is shown in both an example with independent discrete variables and with
continuous correlated variables.

  * •

The same risk measure should be used for the portfolio and in the capital
allocation.

  * •

Simulation may be used in estimating the expected RORAC over the commitment
period of an asset.

  * •

Markov chain Monte Carlo may be used in the estimation of the capital
allocation.

## References

  * [1] Acerbi, C. and Tasche, D. Expected shortfall: a natural coherent alternative to value at risk, Economic notes, 31(2) 376-388, 2002.
  * [2] Artzner, P., Delbaen, F., Eber, J.-M. Heath, D., Thinking coherently, Risk 10 (11) 1997.
  * [3] Artzner, P., Delbaen, F., Eber, J.-M. Heath, D., Coherently measures of risk, Math. Fin. 9 (3) 1999.
  * [4] Denault, M., Coherent allocation of risk capital, Journal of risk 4(1), 1-34, 2001
  * [5] Fischer, T., Risk capital allocation by coherent risk measures based on one-side moments, Insurance: Mathematics and Economics, 32(1), 135-146, 2003
  * [6] Kalkbrener, M., Lotter, H., Overbeck, L. Sensible and efficient capital allocation for credit portfolios, RISK 17(1):S19:24, 2004.
  * [7] Litterman, R., Hot spots and hedges, The Journal of portfolio management, 22, 52-75, 1996.
  * [8] Meyn, S. O and Tweedie, R. L., Markov chains and Stochastic stability, Springer-Verlag, 1993.
  * [9] Tasche, D., Risk contribution and performance measurements, working paper, Technische Universitet Munchen, 1999.
  * [10] Tasche, D., Expected shortfall and beyond, Journal of banking and finance 26(7) 1519-1533, 2002.
  * [11] Tasche, D., Capital Allocation with CreditRisk+, In: CreditRisk+ in the Banking Industry, V.M. Gundlach and F.B. Lehrbass (eds.), pp 25-44. Springer, 2004,
  * [12] Tasche, D., Allocating portfolio economic capital to sub-portfolio, In: Economical capital: a practitioners guide A. Dev (ed.), Risk books, pp. 273-302, 2004.
  * [13] Tasche, D., Measuring sectoral diversification in an asymptotic multi-factor framework, Journal of Credit risk 2(3) 33-55 (2006)
  * [14] Tasche, D. Euler allocation: Theory and practice, Working paper, 2007.
  * [15] Urban, M, Dittrich, J., Kluppelberg, C., and Stolting, R., Allocation of risk capital to insurance portfolios, Blatter der DGVFM, 26, 389-406, 2004