1. No category

Effektiv risikoberegning af aktieporteføljer ved
hjælp af principal-aktiver
Nikola Schou
3. november 2014
Resumé
Der indføres et nyt begreb kaldet principal-aktiver for at muliggøre mere effektive risikoberegninger af aktiemarkeder. Begrebet bygger
på principal-komponent analyse og der argumenteres for at principalaktiverne udgør en optimal måde at repræsentere positioner i et aktiemarked - optimal i forhold til at udregne forskellige risikomål mest
effektivt og optimal i forhold til at analysere markedets korrelationer,
korrektioner, afkast, regimer med mere. For at underbygge teorien gennemføres en række empiriske undersøgelser af amerikanske aktiemarkeder og det vises at beregningstiden for eksempelvis varians, value-at-risk
og expected shortfall kan reduceres med op til 99%.
Vejleder: Kristian Sørensen
1
Indhold
1 Introduktion
4
1.1
Motivation og emnekreds . . . . . . . . . . . . . . . . . . . . .
4
1.2
Problemformulering . . . . . . . . . . . . . . . . . . . . . . . .
7
1.3
Opgavens disposition . . . . . . . . . . . . . . . . . . . . . . .
9
1.4
Tak til . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2 Matematisk baggrund
10
3 Markedet
13
3.1
Kildedata . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3.2
Datasæt og perioder . . . . . . . . . . . . . . . . . . . . . . .
14
4 Principal komponent analyse
14
4.1
Introduktion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
4.2
Matematisk udledning . . . . . . . . . . . . . . . . . . . . . .
16
4.3
Udregning af efficient rand vha. principal-aktiver . . . . . . . .
20
4.4
Empiriske resultater . . . . . . . . . . . . . . . . . . . . . . .
22
5 Effektive beregningsmetoder
24
6 Approximation af varians
24
6.1
Bestemmelse af reduktionsfaktor . . . . . . . . . . . . . . . . .
28
6.2
Empiriske resultater . . . . . . . . . . . . . . . . . . . . . . .
29
6.2.1
Historisk std.afv. for en portefølje . . . . . . . . . . . .
29
6.2.2
Valg af porteføljer . . . . . . . . . . . . . . . . . . . . .
32
6.2.3
Resultat af empiriske undersøgelser . . . . . . . . . . .
33
7 Definition af value-at-risk og expected shortfall
34
7.1
Historisk VaR og ESF . . . . . . . . . . . . . . . . . . . . . .
37
7.2
Empiriske resultater . . . . . . . . . . . . . . . . . . . . . . .
37
8 Approximation af value-at-risk
8.1
40
Bestemmelse af VaR i halen . . . . . . . . . . . . . . . . . . .
2
43
8.2
Empiriske resultater . . . . . . . . . . . . . . . . . . . . . . .
9 Approximation af expected shortfall
9.1
Empiriske resultater . . . . . . . . . . . . . . . . . . . . . . .
44
47
47
10 Bestemmelse af cut-off faktor vha. af random matrix teori (RMT)
52
11 Forecasting
57
11.1 Forecasting i Danske Capital - et resume . . . . . . . . . . . .
57
11.2 Forecasting med principal-aktiver . . . . . . . . . . . . . . . .
59
12 Konklusioner og resultater
62
13 Perspektivering
62
Litteratur
64
Appendices
65
A Programkode
65
3
1 Introduktion
1.1 Motivation og emnekreds
Prisudviklingen for Citigroup Inc. i perioden 1977-2014 ligner en bjergformation skabt med dynamit hvilket kan ses på figuren herunder. Fra starten af oktober 2007 sker et dramatisk fald i værdiansættelsen fra godt 500$ i sommeren
2007 til 10$ i marts 2009. Det svarer til en reduktion på ca. 98% af værdien på
blot 18 måneder, hvilket igen svarer til en halvering af værdien ca. hver 3. måned.
I juli måned 2007 var der stadig gang i festen hos Citigroup, hvilket følgende
udtalelse fra Citigroups administrerende direktør illustrerer:
When the music stops, in terms of liquidity, things will be complicated. But as long as the music is playing, you’ve got to get up
and dance. We’re still dancing. Citigroup Inc. CEO Chuck Prince,
July 2007
På dette tidspunkt vidste hverken ledelsen i Citigroup eller resten af verden at
der var noget så fundamentalt galt med virksomhedens dispositioner, at man
blot 3 måneder senere skulle være vidne til starten på et af de bratteste fald
fra de øverste tinder til fallitens rand i Wall Streets historie.
Før den finansielle nedsmeltning i 2007-2009 var Citigroup en af de største
banker i USA med 357.000 ansatte1 . Da Citigroup i perioden op til oktober
2007 havde status som ledende børsnoteret selskab i USA, må man antage at
adskillige kloge hoveder i og udenfor virksomheden har lavet daglige analyser
af Citigroup med henblik på at værdiansætte virksomhedens aktier, herunder
at vurdere risikoen i de aktiver, banken lå inde med. Hvordan kan man tage
1
Se http://en.wikipedia.org/wiki/Citigroup
4
så gruelig fejl, at alle - inklusiv den øverste administrerende direktør - i sidste
ende fremstår som glade amatører, der blot danser så længe musikken spiller
uden en klar idé om, hvornår den stopper eller hvor mange stole der vil være
tilbage at sidde på?
Vurdering af risiko er en videnskab men også - som eksemplet med Citigroup
illustrerer - konsensuspræget gætværk på kanten til kollektivt bedrag. Ikke
desto mindre er risikovurdering en grundpille i alle moderne samfund, hvilket
er motivationen for at denne opgave først og fremmest skal handle om risiko.
Ud af stor risiko kommer stor fortjeneste men også store tab. I et samfundsperspektiv er store fortjenester til enkelte individers fordel, risikotagerne, mens
store tab potentielt kan føre til andre endnu større tab og i sidste ende destabilisere økonomien med negative konsekvenser for langt flere personer, herunder
personer som ingen involvering i øvrigt har i den risiko, der blev taget, ej heller med udsigt til deling af den potentielle fortjeneste. Denne asymmetri giver
anledning til det velkendte begreb moral hazards, hvor risikotagerne fristes til
at tage en for høj risiko, ganske simpelt fordi de ved, at en potentiel optur er
til egen vinding mens en potentiel nedtur deles med offentligheden.
Af disse årsager er risikostyring et vigtigt nationalt og internationalt samfundsanliggende, som har udmøntet sig i en række finansmyndigheder, tilsyn og regulativer. Konsekvensen heraf er, at der stilles stadig større krav til
risikostyringen i alle større virksomheder, herunder i særlig grad finansielle
virksomheder såsom banker, investeringsforeninger og pensionskasser. Sådanne virksomheder er i dag nødt til at besidde folk med særlige kompetencer og
specialiseret software, der tilsammen giver virksomheden en effektiv forståelse
og kvantificering af den risiko, de står med, herunder skal de beherske effektive metoder til at lave fremskrivninger af volatilitet, afkast og andre størrelser
med relevans for risikostyring og/eller porteføljestyring.
Forskellige virksomheder er udsat for forskellige former for risici. Eksempelvis
er produktionsvirksomheder især udsat for operationel risiko mens finansielle
virksomheder især er udsat for markedsrisiko, equity risk, kreditrisiko et.al.
Forskellige former for risiko kræver forskellige modeller, både konceptuelt og
beregningsmæssigt.
For at analysere risiko er der i de seneste 5-10 år opstået interesse for et begreb
kaldet risikofaktorer, se [13] og [8]. En risiko-faktor er basalt set en makrostørrelse i markedet såsom markedsværdi, momentum, volumen og volatilitet.
På Figur 1 ses eksempler på risiko-faktorer samt korrelationer imellem dem.
Disse makro-størrelser kan kvantificeres og der kan laves lineær regression
af de enkelte aktiver i markedet mod disse makro-størrelser. Hermed ender
man med at kunne opskrive ethvert aktiv eller portefølje i markedet som en
linearkombination af risiko-faktorerne. Ifølge forfatterne til [8] får man hermed
en forbedret model til vurdering af risiko og allokering af efficiente porteføljer.
5
Figur 1: Tabel kopieret fra [8]. Viser eksempler på risikofaktorer samt
deres interne korrelationer.
I denne opgave tager vi skridtet videre ved at indføre begrebet principalaktiver. Principal-aktiver kan opfattes som en slags risikofaktorer der dog adskiller sig fra andre faktor-modeller på følgende punkter:
• Principal-aktiverne fremkommer via en matematisk velfunderet model
og ikke ud fra et gæt om hvad de vigtige makro-økonomiske størrelser
er.
• Principal-aktiverne udgør en såkaldt fuldstændig basis til at beskrive
positioner (dvs. porteføljer) i markedet. Hermed menes at alle positioner kan repræsenteres eksakt og entydigt som en linear-kombination af
principal-aktiverne. Samtidig har principal-aktiverne en orden, hvor det
første principal-aktiv er det vigtigste, det andet er det næstvigtigste og
så fremdeles. Det betyder, at man kan lave graduerede repræsentationer
af en given position, der spænder lige fra anvendelse af blot det første
principal-aktiv - som giver den højeste beregnings-effektivitet men den
dårligste approximation - til anvendelse af alle principal-aktiver, der giver den laveste beregnings-effektivitet men til gengæld opretholder fuld
præcision.
• Principal-aktiverne er ukorrelerede. Dette er i modsætning til normale
risiko-faktorer der har korrelationer, jævnfør Figur 1.
• Principal-aktiverne gør det muligt at approximere forskellige risikomål
såsom varians, value-at-risk og expected shortfall på en overraskende
effektiv måde.
Disse egenskaber gør det yderst attraktivt at beskæftige sig med principalaktiver og dermed er der redegjort for, at dette begreb udgør fundamentet i
denne opgave.
6
Et berømt udsagn lyder:
Time is money. - Benjamin Franklin
Sagt med andre ord:
Time is money says the proverb, but turn it around and you get
a precious truth. Money is time. - George Gissing
Time is money. Wasted time means wasted money means trouble.
- Shirley Temple
Tid og penge er yin og yan indenfor økonomisk virke. Og samtidig er penge
og dermed tid altid den vigtigste knappe ressource.
Risikovurderinger for finansielle virksomheder og især risikovurderinger, der
skal produceres til finansmyndighederne, handler i sidste ende om at lave nogle
udregninger, der ender med kvantitative udsagn formuleret med tal eller grafer. Fælles for alle disse risikoberegninger er, at de næsten altid baserer sig på
tidsserier for historisk udvikling af forskellige størrelser. Dette giver meget ofte
anledning til ekstreme datamængder, som skal indsamles, lagres og behandles. Og eftersom risikovurderingerne skal laves hver dag eller hver uge, ender
mange risikoberegninger med at være meget tidskrævende. Af denne grund
er effektive beregningsmetoder - eller algoritmer om man vil - værdifulde for
store virksomheder. En effektiv approximation, der giver en reduktion i en given risikoberegning med 99% kan eksempelvis gøre forskellen på, om man kan
producere risikoberegninger i real-tid eller ej. Eller på om ledelsen kan nå at
reagere i rette tid på en given hændelse. Derfor er effektive beregningsmetoder
for risiko guld værd for store finansielle virksomheder.
For at skabe en overkommelig ramme for opgaven, vil vi begrænse os til risiko
på aktiemarkeder og hermed har vi redegjort for det nøjagtige ordvalg i opgavens titel: Effektiv risikoberegning af aktieporteføljer ved hjælp af principalaktiver
1.2 Problemformulering
Hvordan analyseres et aktiemarked mest effektivt?
Hvordan ser vi klarest, hvornår der sker korrektioner i markedet?
Hvordan slipper vi ud over den vilkårlighed der ligger i at udvælge markedsindekser eller risiko-faktorer baseret på makro-økonomiske størrelser?
7
Kan man forestille sig et analyse-apparat der reducerer kompleksitet og beregningstid for risikoberegninger radikalt?
Denne opgave vil forsøge at give svar på disse spørgsmål. Som nævnt flere
gange begrænser vi os til aktiemarkeder. De resultater og metoder der udvikles
kan i en vis udstrækning anvendes på alle aktiver, der handles på et marked
med tilstrækkelig høj likviditet og efficiens, men i visse sammenhænge gøres
der antagelser, som holder bedst for diversificerede aktieporteføljer.
I traditionelle faktor-modeller laves der som allerede nævnt lineær regression
imod forudbestemte makro-størrelser. Vi vil i stedet tage udgangspunkt i en
stringent matematisk metode, som ligger i grænselandet mellem statistik og
lineær algebra (eller teori om mange-dimensionelle vektorrum om man vil).
Metoden hedder principal-komponent analyse og er vidt og bredt beskrevet
og endvidere anvendt indenfor finansiering igennem mange år, se fx. [1] afsnit
4.5 "Interpreting principal components - Stock Market Prices"eller [2] afsnit
3.4.4 eller [3] side 208.
Denne metode giver den fundamentale indsigt, at der i et marked med n
aktiver entydigt må findes n kombinationer af aktiver som har den egenskab
at de (i) er fuldstændig ukorrelerede, (ii) hver især maksimerer variansen og
(iii) udspænder markedet. Disse tre egenskaber kræver naturligvis en del mere
forklaring, og det vil blive givet i afsnit 4. Selv uden yderligere forklaring vil
nysgerrigheden hos mange læsere nok være vakt i forhold til hvordan disse
kombinationer af aktiver findes, hvordan de ser ud og hvad de kan bruges
til. I opgaven vil vi kalde dem for principal-aktiver. Hermed understreger vi
tilhørsforholdet til principal-komponent analysen samt det faktum, at vi kan
betragte dem som et nyt aktiv på linie med både de oprindelige aktiver og
alle andre mulige porteføljer i markedet.
Den overordnede tese for opgaven er, at principal-aktiverne giver et optimalt
eller tæt-på-optimalt udgangspunkt for at analysere risiko i markedet, herunder at lave effektive approximationer og beregningsmodeller for diversificerede
aktieporteføljer.
Som centralt resultat skal nævnes approximationer og beregningsmetoder for
historisk standard-afvigelse, historisk value-at-risk og historisk expected shortfall, med reduktioner på op til 99% af beregningstiden. Opgaven vil både udlede teorien bag og gennemføre en empirisk efterprøvning af approximationernes
validitet.
Herudover vil opgaven give en introduktion til principal-komponent analyse
uden den traditionelle brug af Lagrange multipliers. Derimod tager vi udgangspunkt i et centralt resultat fra lineær algebra, den såkaldte spektralsætning.2
Alle relevante egenskaber ved principal-aktiverne udledes fra denne læresæt2
Spektralsætningen for n-dimensionelle vektorrum følger af egenskaber ved determinanter og n’te-grads polynomier. Et kortfattet bevis kan findes her: http://tinyurl.com/ljvvypq
8
ning med relativt få matematiske forudsætninger. Denne fremgangsmåde er
valgt istedet for Lagrange multipliers, dels fordi bevisførelsen efter min mening er lettere og mere elegant, dels for at tilbyde et alternativ til traditionelle
lærebøger i finansiering eller principal komponent-analyser. Hvis læseren ønsker at se en traditionel fremstilling via Lagrange multipliers henvises til [1],
afsnit 1.1.
For yderligere at underbygge de approximationer der gøres, inddrages udvalgte resultater fra random matrix-teori. Ved hjælp af den såkaldte MarcenkoPastur fordeling argumenteres der for, hvor mange principal-aktiver der har
statistisk signifikans.
Slutteligt vil opgaven i forbindelse med forecasting introducere nogle nye værktøjer til analyse af markedskorrektioner baseret på principal-aktiverne og afledte begreber. Vi vil argumentere for, at principal-aktiverne og relaterede
teknikker giver en god forståelse for markedsdynamikkerne. Særligt vil det
blive forklaret hvorledes korrektioner i markedet kan afkodes forholdsvis klart,
hvilket i sidste ende gør det muligt at lave bedre modeller for forecast af risiko.
1.3 Opgavens disposition
• Den nødvendige og tilstrækkelige matematiske baggrund for opgaven
gennemgås i afsnit 2.
• Markedet og kildedata beskrives i afsnit 3. Dvs. der redegøres for hvilket
marked vi befinder os på, hvilke datasæt der arbejdes med, hvor data
hentes fra og hvordan det er blevet bearbejdet.
• Afsnit 4 giver en introduktion til principal-komponent analysen med
en udledning af relevante resultater. Samtidig introduceres principalaktiver som begreb.
• Afsnit 5 formaliserer idéen om effektive beregningsmetoder. Herunder
gives en simpel måde at måle hvor stor reduktion i beregningstid en
given approximation giver.
• I afsnit 6 udledes en effektiv approximation af historisk varians/std.afvigelse
med udgangspunkt i principal-aktiverne. Herunder vises det at man kan
opnå en reduktion i beregningstiden af historisk varians/std.afvigelse på
90-99%.
• Value-at-risk og expected shortfall defineres som risikomål i afsnit 7.
Samtidig beskrives hvorledes de kan bestemmes historisk ud fra en given
tidsserie.
9
• I afsnit 8 udledes en effektiv approximation af historisk value-at-risk og
expected shortfall. Igen vises det at beregningstiden kan reduceres med
90-99%.
• I afsnit 9 gentages øvelsen fra afsnit 8, nu blot for expected shortfall.
• Afsnit 10 inddrager nogle centrale resultater fra random matrix theory
til bestemmelse af signifikante og ikke-signifikante principal-aktiver.
• Afsnit 11 vil med udgangspunkt i nogle modeller fra Danske Capital
introducere forecasting som emne og forklare hvorledes principal-aktiver
kan anvendes indenfor dette felt.
• Afsnit 12 vil opridse de vigtigste konklusioner mens afsnit og 13 giver
en yderligere perspektivering af principal-aktiver og deres anvendelsesmuligheder.
Der er til denne opgave anvendt en række artikler, som findes på internettet. For at gøre det nemt for læseren kan de hentes samlet i denne
zip-fil:
http://nikolaschou.dk/principalaktiver/artikler.zip
1.4 Tak til
Tak til Kristian Sørensen, CBS, for omhyggeligt vejlederskab, til Peter Dalgaard, CBS, for at give vejledning vedr. statistisk teori, til Rasmus Majborn
og Jacob Hannover fra Danske Capital for at stille op til interview og udlåne
materialer, til Erik Gadeberg og Anders Haaning fra Jyske Bank for at perspektivere emnet og til Scott Ehlers og Else Braathen fra SimCorp for at stille
op til et interview tidligt i opgavens forløb.
2 Matematisk baggrund
Til denne opgave anvendes følgende notationer, konventioner og relevante regneregler indenfor matrixregning og statistik.
• Vektorer angives med små bogstaver i fed skrift, fx. w, hvor elementerne
i vektoren angives w1 , w2 , . . . , wn .
• Vektorer kan indekseres hvis der er m af slagsen, fx w1 , w2 , . . . wm . I
dette tilfælde angives elementerne af hver enkelt vektor med et ekstra
index, fx w1 = (w11 , w12 , . . . , w1n )
• Matricer angives med store bogstaver i fed skrift, fx. C, hvor de enkelte
elementer i C angives som Cij .
10
• Matricer og vektorer kan som bekendt transponeres og dette angives
med et 0-symbol, dvs.
a b
A= c d
a
⇐⇒ A0 = b
c
d
og

w=
w1 w2

w1
 w2 

. . . wn ⇐⇒ w0 = 
 ... 
wn
• Vektorer er som udgangspunkt rækkevektorer. Ved at transponere en
vektor w får man en søjlevektor w0.
• Matrix-transponering opfylder følgende vigtige regneregler: (A0 )0 = A
og (AB)0 = B 0 A0
• Matrix-multiplikation - herunder multiplikation mellem vektorer og matricer - er associativ, dvs. A(BC) = (AB)C. Det understreges at matrixmultiplikation ikke er kommutativ, dvs. AB er generelt set ikke det
samme som BA.
• Sporet af en matrix (på engelsk trace) defineres som summen af diagonalelementerne:
X
Tr(A) =
Aii
i
˙ herunder at Tr(AB) =
• Der gælder forskellige praktiske regneregler for Tr(),
Tr(BA).
• For en vektor w betegner diag(w) den matrix som har w1 , . . . , wn langs
diagonalen og 0 for resten. For en matrix C betegner diag(C) den vektor
som består af diagonalelementerne (C11 , . . . , Cnn ). Man kan sige at diag()˙
afbilder frem og tilbage mellem vektorer og matricer.
• I n defineres som enhedsmatricen i Rn , dvs. I n = diag(1, 1, . . . , 1).
• Stokastiske variable angives med ’krøllede’ bogstaver, fx A, B, C, P, Q
osv. Flere stokastiske variable indekseres som A1 , . . . , An og dette kan
knyttes sammen i en vektor angivet med en overstregning, fx A =
(A1 , . . . , An ).
11
• For en stokastisk variabel A og en fraktil α mellem 0 og 1 gælder det,
at E(A) og µA angiver middelværdien„ Var(A) angiver variansen, σ(A)
angiver std.afvigelsen (spredningen), ESFα (A) angiver expected shortfall
ved α og VaRα (A) angiver value-at-risk ved α. Bemærk at der gives en
detaljeret introduktion og definition af value-at-risk og expected shortfall
i afsnit 7.
• I klassisk vektorregning definerer man skalar-produktet mellem a og b
(også kaldet det indre produkt) som summen af produktet af elementerne. I denne opgave anvendes udelukkende matrix-multiplikation for at
holde en simpel notation, dvs. hvis vi ønsker at finde skalar-produktet
mellem to række-vektorer gøres det via matrix-operationen: ab0 som
alternativt skrives a · b0 .
• Hvis man multiplicerer a med b0 får man som nævnt skalar-produktet.
Hvis man derimod multiplicerer b0 med a får man en n × n matrix,
hvilkes ses her:
( a
 
 1
b 1 b 1 a1

 b 2 a1
b
2
 .
b0 a ≡ b0 · a = . 
 .. 
 ..
bn
b n a1
a2
b1 a2
b2 a2
..
.
...
...
...
...
an )

b1 an
b2 an 
.. 
. 
b n a2
...
b n an
• Længden af en vektor defineres som
v
!
u n
u X
a2i
kak = t
i=1
• Vinklen θ mellem to vektorer a og b findes fra formlen:
cos(θ) = kakkbka · b0
• En enhedsvektor er en vektor som har længden 1.
• De enhedsvektorer som

1
 0


i1 =  0
 ..
 .
0
peger langs akserne betegnes

 

0
0

 1 
 0

 


 0 
 .
 , i2 =   , . . . , in =  ..

 .. 

 . 
 0

0
1







• To vektorer a og b er ortogonale eller vinkelrette hvis a · b0 = 0 hvilket
let ses ved θ = cos−1 (0) = π/2 = 90◦ .
12
3 Markedet
Vi betragter et marked bestående af n aktier og lader A1 , . . . , An betegne
stokastiske variable, som angiver den daglige relative prisændring for de n
aktier.
Vi antager at den gennemsnitlige prisændring henover tid (den langsigtede
vækst) er meget mindre end de daglige udsving eller vi kan fratrække den
langsigtede vækst fra de daglige udsving, dvs. vi antager E(Ai ) ≡ µAi = 0 for
alle i.
I dette marked betragter vi en portefølje P karakteriseret ved en vægtet sammensætning af de n aktier, hvor vægtene betegnes w = (w1 , . . . , wn ). Dvs. vi
kan tænke på P som en stokastisk variabel defineret som følger:
P=
n
X
0
wi Ai = wA hvor
X
wi = 1
(1)
i
i=1
For en nemheds skyld forestiller vi os, at alle aktier kan
Phandles i enheder, der
hver koster 1kr hvormed også P koster 1kr, eftersom i wi = 1.
Det bemærkes at (1) implicit antager, at porteføljen løbende rebalanceres.
Hvis fx den absolutte pris på Ai stiger mens alle øvrige er uændrede, da vil
man være nødt til at sælge ud af Ai og købe op i de øvrige for at fastholde en
given sammensætning af vægte.
3.1 Kildedata
Som datakilde er anvendt Yahoo Finance, et gratis web-API udstillet af Yahoo.
Data hentes og analyseres vha. Matlab.
I princippet er det simpelt at hente data fra Matlab via funktionen buildUniverse()3 . Problemet er blot at de simple metoder ikke virker idet at priserne
ikke er korrigeret i forhold til dividender og aktiesplit. Det betyder at man reelt
ikke kan bruge disse data til at finde daglige afkast. Istedet er man nødt til at
bruge funktionen fetch(), som giver mulighed for at hente såkaldte ’Adjusted
Close’ priser.
Næste udfordring er at finde de sæt af aktiver, man ønsker at analysere. På
den ene side vil man gerne dække så mange aktiver som muligt. Samtidig vil
man gerne dække så lang en periode som muligt. Det leder en til spørgsmålet
om, hvorledes man kombinerer lange tidsserier med korte, dvs. hvorledes laves
analysen når nogle aktiver har eksisteret siden 60’erne, andre først er kommet
til senere. Til denne opgave er valgt følgende konventioner for en analyse, der
3
Dokumentation af denne funktion kan læses her:
http://www.mathworks.se/help/datafeed/yahoo.builduniverse.html?refresh=true
13
løber over en given periode: (i) Når et aktiv introduceres på børsen med en
pris P på en dato D efter periodens start, vil prisen i intervallet [start, D] blive
sat til P. (ii) Når et aktiv fjernes fra børsen sker det samme, dvs. den sidste
pris får lov at fortsætte uændret.
Sidste udfordring bunder i, at forskellige aktiver handles på forskellige handelsdage. Det kan enten skyldes at aktiverne handles på forskellige børser,
eller det kan skyldes at Yahoo har ufuldstændige data - hvilket desværre er et
faktum4 . Det gør det en smule kompliceret at opstille en samlet tidsserie, da
det kræver en form for sammenfletning. Tænk blot på en tidsserie A som har
priser for 2, 3 og 6. januar (i et givet år) og en anden tidsserie B som har priser
for 2, 4, og 5. januar. For denne analyse er det valgt at tage udgangspunkt i
de handelsdage, der er registreret for IBM. Alle andre tidsserier er flettet ind
således at priserne listes på IBM’s handelsdage vha. en slags best-fit algoritme.
Koden kan ses i appendix og er markeret med teksten ’MERGING TRADING
DAYS’.
3.2 Datasæt og perioder
Alle empiriske undersøgelser er baseret på et eller flere af de datasæt der
vises i Tabel 1. For alle tre datasæt har det været nødvendigt at rense ud i
de anvendte aktiver af forskellige årsager. Den primære årsag har været, at
Yahoo Finance ikke kunne levere valide data for en del af aktiverne. Således
var det for datasættet NYSE77 nødvendigt at fjerne ca. 177 aktiver for at få
et rent datasæt. Metoden til at rense ud i datasættene har været, at finde
de aktiver med særlig usædvanlige spring i prisudviklingen (fx. en stigning
på 5000% på en dag). Dernæst har det været en manuel opgave at gå dem
igennem og be- eller afkræfte, hvorvidt de leverede data er korrekte eller ej.
Vurderingen af om de er korrekte er taget ved at sammenligne priskurven med
den officielle priskurve fra de grafer man kan finde online via Yahoo Finance
web-grænsefladen.
4 Principal komponent analyse
4.1 Introduktion
Principal komponent-analyser (forkortet PCA) som metode blev opfundet af
Karl Pearson i 1901 og har siden fundet anvendelser indenfor mange viden4
For en diskussion af fejl og mangler i Yahoo Finance data se her:
http://quant.stackexchange.com/questions/942/any-known-bugs-with-yahoo-financeadjusted-close-data
og her:
http://quant.stackexchange.com/questions/7216/daily-returns-using-adjusted-close
14
Tabel 1: Oversigt over datasæt anvendt til de empiriske undersøgelser
ID
Beskrivelse
SP2000 Indeholder de aktier fra S&P500
indekset som har været i markedet siden år 2000
SP1984 Indeholder de aktier fra S&P500
indekset som har været i markedet siden år 1984
NYSE77 Indeholder de aktier fra New York
Stock exchange som har været i
markedet siden år 1977
# aktiver
n
463
Periode
2000-2014
#dage
T
3521
98
1984-2014
7565
1594
1977-2014
13090
skabelige discipliner, hvor metoden har fået forskellige navne, fx egenværdi
dekomposition indenfor lineær algebra, spektralanalyse indenfor akustik og fysik med mere, jf. [14].
For at få den rigtige intuition omkring PCA indenfor rammerne af denne opgave, kan det være en god ide at starte med at diskutere markedsindekser [12].
Markedsindekser dannes som bekendt ved at lave en fiktiv portefølje bestående af bestemte aktier med nogle bestemte vægte, hvormed man løbende kan
udregne en indeks-pris ud fra vægtene og de underliggende aktiepriser. Hermed får man en simpel én-dimensionel størrelse som er nemmere at forholde
sig til end markedet som helhed.
Der findes en lang række forskellige markedsindekser med forskellige kriterier
for udvælgelse af aktier og forskellige metoder til bestemmelse af vægtene.
Nogle indekser forsøger at samle aktier fra alle de største virksomheder på
tværs af hele markedet, globalt eller nationalt. Fx samler S&P Global 100indekset 100 multinationale aktier. Andre indekser samler aktier indenfor et
bestemt land, en bestemt sektor, en bestemt børs osv. osv.
Med hensyn til valg af vægte er der primært to metoder i anvendelse. For
prisvægtede indekser vælges vægten ud fra prisen. Hvis en aktie er dobbelt
så dyr som en anden aktie, vejer den også dobbelt så meget i indekset. For
kapitalvægtede indekser er det virksomhedens samlede markedsværdi, der bestemmer vægten. Hvis én aktie har dobbelt så stor markedsværdi som en anden
aktie, da er vægten også dobbelt så stor [14].
Markedsindekser bruges til mange forskellige slags formål, men nok primært
til markedsanalyser. Hvis et indeks stiger ses det som et udtryk for en generel
stigning i markedet. Hvis et indeks har kraftige udsving ses det som udtryk
for, at markedet generelt har høj volatilitet. Hvis udviklingsmønstrene i indekset ændrer sig, ses det måske som et udtryk for et regimeskifte i markedet.
Indekserne anses altså for at være gode aggregerede størrelser, hvorudfra vi
15
kan aflæse og forstå markedet.
Man kan i denne sammenhæng undre sig lidt over den vilkårlighed der er i
konstruktionen af indekserne. Og det kan virke uoverskueligt eller overflødigt
at der er så mange indekser at vælge imellem. Principal komponent-analysen
kan motiveres ud fra et ønske om at frembringe et mindre antal indekser som
kan bygges på et teoretisk grundlag, både matematisk og finansielt. Med dette
som indgangsbøn udleder vi nu principal komponent analysens grundsætninger og en række afledte værktøjer.
4.2 Matematisk udledning
Vi starter med at minde om at vi betragter et marked med n aktier betegnet
A1 , . . . , An og herudfra kan konstrueres porteføljer som en vægtet kombination
af aktierne, hvor vægtene betegnes med en vektor w. Vi minder også om at
vi i vores model tænker på både aktier og porteføljer som stokastiske variable
der angiver det daglige relative afkast, og at en portefølje kan skrives som
0
matrixprodukt: P = wA .
Vi starter nu med at udregne variansen for P.
0
0 0
2
Var(P) = E((P) ) = E wA wA
0
0
= E(wA Aw0 ) = wE(A A)w0
Her er det udnyttet at middelværdien E(·) som funktion er lineær. Nu indfører
man kovarians-matricen C defineret ved:
0
C ≡ Cov(A) ≡ E(A A) eller alternativt Cij = E(Ai Aj ) for alle i, j
og vi kan opsummere følgende vigtige sammenhæng:
Var(P) = wCw0 hvor Cij = E(Ai Aj ) for alle i, j
(2)
Det ses at kovarians-matricen er symmetrisk, dvs. Cij = Cji . Vi vil nu udnytte
det såkaldte spektral teorem, der er en hovedsætning indenfor lineær algebra og
matrixregning.5 Spektral-teoremet siger, at enhver symmetrisk matrix C
kan diagonaliseres, dvs. vi kan omskrive C som et produkt af tre matricer
som følger:
5
Spektral teoremet kan også formuleres mere generelt og indenfor mere avancerede grene
af matematikken, men til denne opgave er vi kun interesseret i en mere snæver anvendelse.
16
C = QΛQ−1
(3)
hvor Q er en ortonormal matrix (se nedenfor) og hvor Λ er en diagonalmatrix

λ1 0 0 0
 0 λ2 0 0

Λ=
 0 0 ... 0
0 0 0 λn





med faldende værdier langs diagonalen. Dvs. λ1 er den største værdi, λ2 er den
næststørste værdi osv. Der gælder endvidere at denne opsplitning eller dekomposition er entydig. Diagonalværdierne i Λ kaldes traditionelt for egenværdier
mens søjlevektorerne i Q kaldes for egenvektorer eller principal komponenter.
6
NB: Bemærk at ’egenværdi’ og ’egenvektor’ på engelsk kaldes for ’eigenvalue’
og ’eigenvector’ respektivt. Da visse figurer i denne opgave har engelske
figur-tekster, vil disse betegnelser indimellem blive brugt.
En ortonormal matrix består pr. definition af n enhedsvektorer, der alle er
ortogonale (dvs. vinkelrette) på hinanden, jævnfør afsnit 2.
Lad os i det følgende bruge lidt tid på at opbygge lidt intuition omkring C,
Q og Λ og deres indbyrdes sammenhænge.
Det er praktisk at benævne søjlerne i Q med q01 , . . . , q0n , dvs. vi kan skrive
således:
 


   
 


 





0
0
0
≡
q1
q2
qn
...
Q≡

 





 
q11 
 
 q21
..
.
.   ..
q   q
1n
2n





 qn1
..
...
.



 q
nn
hvor de krøllede parenteser kun er medtaget for at markere søjlestrukturen.
At Q er ortonormal betyder per definition at:
qi q0j
=
hvis i 6= j
hvis i = j
0
1
6
Dette udsagn er bevidst forsimplet og bygger på en række antagelser, nemlig at
det(C) 6= 0 og at alle diagonalværdierne er forskellige. For typiske finansielle tidsserier
med flere målepunkter end aktiver vil disse antagelser dog holde.
17
 






hvilket alternativt kan skrives:
QQ0 = Q0 Q = I n hvor I n er enhedsmatricen
Bemærk især fra disse ligninger at transponeringen af Q er sin egen inverse,
dvs.
Q−1 = Q0
(4)
Herfra konkluderes at Q er længdebevarende, hvilket demonstreres af følgende
udregning:
kQa0 k2 = (Qa0 )0 (Qa0 ) = a(Q0 Q)a0 = aI n a0 = aa0 = kak2
Ydermere kan det ses at Q er vinkelbevarende, dvs. hvis θ1 betegner vinklen
mellem a og b og θ2 betegner vinklen mellem a og b efter multiplikation med
C, da finder man:
cos(θ2 ) =
aQ0 Qb0
ab0
(Qa0 )0 Qb0
=
=
= cos(θ1 )
kQakkQbk
kakkbk
kakkbk
En matrix der både er længdebevarende og vinkelbevarende ved multiplikation, kender man også som en rotation. Samtidig ser man at Λ ved multiplikation virker som ren skalering. Vi kan altså betragte produktet QΛQ0 som
en samlet operation der først laver en rotation Q0 , dernæst en skalering Λ og
slutteligt en modsatrettet rotation Q.
Hernæst kan vi indse at følgende gælder:
Cq0i = QΛQ0 q0i = QΛii = Qλi i1 = λi q0i
Med andre ord, C virker på qi som en skalering med skaleringsfaktor λi .
Generelt set vil en matrix udvirke både skalering og rotation når den ganges
på en vektor, men disse særlige vektorer udmærker sig altså ved kun at blive
skaleret og ikke roteret af C. I den forstand kan man sige, at der er et særligt
tilhørsforhold mellem C og qi , hvilket er årsagen til at man kalder q1 , . . . , qn
for egenvektorer for C og de tilhørende λ’er for egenværdier. Hvis vi tænker på
vektorerne som komponenter ledes man til betegnelsen principal komponenter,
heraf forkortelsen PCA for principal komponent analyse.
Den første egenvektor q1 og den første egenværdi λ1 har en endnu mere speciel egenskab. Der gælder nemlig at q1 angiver den enhedsvektor, der giver
maksimal varians efter formel (2). Det kan ses ved at betragte en vilkårlig
enhedsvektor e og finde variansen for den tilhørende portefølje:
18
Figur 2: Illustrerer enhedsvektorerne for n = 3 og ideen om at én bestemt
vektor markeret med tallet 1 udmærker sig ved at maksimere variansen.
Var(Pe ) = eCe0 = eQΛQ0 e0 = (eQ)Λ(eQ)0 = f Λf 0 =
X
fi2 λi
i
≤ λ1
X
fi2 = λ1 kf k2 = λ1 = Var(Pq1 )
i
I udregningen har vi defineret f ≡ eQ og udnyttet at Q er længdebevarende,
dvs. kf k = 1 fordi e har længden 1.
Vi har altså vist at uanset hvordan vi vælger enhedsvektoren e vil det gælde
at Var(Pe ) ≤ Var(Pq1 ). Ideen er illustreret på figur 2.
Ved at anvende en lignende teknik kan man også indse, at blandt alle enhedsvektorer som er vinkelret på q1 , da er q2 den som giver størst varians og mere
generelt, blandt alle enhedsvektorer som er vinkelrette på q1 , q2 , . . . , qk , da er
qk+1 den som giver størst varians.
Lad nu Qi betegne den portefølje som har qi som vægte, dvs.
Qi =
n
X
0
qij Aj eller Qi = qi · A eller Q = Q0 · A
j=1
Da finder man variansen:
Var(Qi ) = qi Cq0i = qi (λi qi0 ) = λi
og kovariansen for i 6= j:
0
Cov(Qi , Qj ) = E(Qi · Qj ) = E qi A Aq0j
0 = qi E A A q0j = qi Cq0j = qi λj q0j = λj (qi q0j ) = 0
19
(5)
(6)
Det vil altså sige at Q1 , . . . , Qn udgør n porteføljer som er gensidigt ukorrelerede og som har variansen givet ved egenværdierne λ1 , . . . , λn .
Herudover gælder det at enhver portefølje P med vægte w entydigt kan skrives
som en linearkombination af Q1 , . . . , Qn som det vises her:
0
0
0
P = wA = (wQ) Q A
0
=β·Q =
n
X
βi Qi hvor β ≡ w · Q
(7)
i=1
hvilket også kan skrives som:
P = β1 Q1 + · · · + βn Qn hvor βi ≡ w · q0i
I denne ligning fortolkes βi som projektionen af w på den i’te egenvektor qi .
Vi ser altså at Q1 , . . . , Qn udspænder markedet og ved at operere med βi ’erne
istedet for wi ’erne får vi åbenbart en alternativ måde at repræsentere og analysere porteføljer i markedet. Fx kan vi finde variansen således:
Var(P) = Var(β1 Q1 + · · · + βn Qn ) = β12 Var(Q1 ) + · · · + βn2 Var(Qn )
Her er det udnyttet at alle krydsledene går ud eftersom Q1 , . . . , Qn er ukorrelerede. Med andre ord:
σP2 = Var(P) = β12 λ1 + · · · + βn2 λn hvor βi ≡ w · q0i
(8)
Det er her udnyttet at Qi ’erne er ukorrelerede hvorved alle ’krydsleddene’ udgår. Samtidig er det udnyttet at variansen af Qi er λi som vist ovenfor. Den
vigtige observation er her, at når en portefølje udtrykkes som en vægtet
sum af de principale komponenter (eller egenvektorerne), da findes
variansen som den vægtede sum af variansen af principal komponenterne. Dette er praktisk og elegant. Vi skal simpelthen tænke på markedet
som om de primære aktiver er Q1 , . . . , Qn og enhver portefølje (herunder de
oprindelige enkelte aktier i markedet) udtrykkes som linearkombinationer af
Qi ’erne. Man kan også sige at Qi ’erne udgør en ny basis for markedet. Denne
ide er illustreret på figur 3.
4.3 Udregning af efficient rand vha. principal-aktiver
Som en anvendelse af de udledte formler for principal-aktiverne vil vi nu se,
hvorledes punkter på den efficiente rand kan udtrykkes ved hjælp af principal20
Figur 3: Illustrerer ideen om at finde en ny basis, som ’passer’ bedre
til det givne problemfelt, her illustreret som en række målepunkter i 2
dimensioner. I den nye basis bestående af y1 og y2 ses målepunkterne at
have en variation der følger akserne. I den gamle basis bestående af x1
og x2 ses målepunkter at ligge på skrå.
aktiver. Vi tager udgangspunkt i den velkendte formel for udledning af markedsporteføljen på den efficiente rand jævnfør [10] ligning (3.20):
0
w
ˆ market
= C −1 (µ0A − r0 )
hvor µA = (µA1 , . . . , µAn ) angiver det gennemsnitlige afkast af aktiverne
A1 , . . . , An , r angiver den risikofri rente som en konstant vektor (dvs. en vektor med samme værdi r på alle pladser) og w
ˆ market angiver den unormaliserede
markedsportefølje. Dvs. den reelle markedsportefølje findes ved normalisering:
0
0
wmarket
=w
ˆ market
1
|w
ˆ market |
hvor | · | betegner summen af elementerne
Ved at bruge (3) og (4) kan vi nu omskrive denne formel på følgende måde:
−1
0
w
ˆ market
= (QΛQ0 )
(µ0A − r0 )
= (Q0 )−1 Λ−1 Q−1 (µ0A − r0 )
= QΛ−1 Q0 (µ0A − r0 )
= QΛ−1 µ0Q − r0Q
hvor vi har indført
µ0Q ≡ Q0 µ0A , dvs. forventet afkast af Q1 , . . . , Qn


|q1 |
 |q2 | 


r0Q ≡ r  ..  fortolket som en normaliseret risikofri rente
 . 
|qn |
21
Ved at gange med Q0 fra venstre på begge sider af lighedstegnet og udnytte
definitionen af β fra (7) når man til konklusionen:
0
βˆmarket = Λ−1 µ0Q − r0Q eller alternativt
βˆmarket,i = λ−1
i (µQi − r|qi |) hvor
0
1
β 0market = βˆmarket
|βˆmarket |
(9)
Vi ser altså, at porteføljer på den efficiente rand kan konstrueres ud fra
principal-aktiverne hvor vægten βi for det i’te principal-aktiv Qi , er proportionel med det forventede afkast af Qi fratrukket den risikofri rente gange prisen
|qi | på Qi og dernæst divideret med variansen λi af Qi .
4.4 Empiriske resultater
For at få en første validering af PCA analysen vil vi finde de første tre
principal-komponenter Q1 , Q2 , Q3 for markedsdata hørende til datasættes SP84,
jævnfør afsnit 3.2. Dernæst vil vi finde prisudviklingen for de tre principalkomponenter. Det gøres ved dag for dag at finde prisen på det indeks som har
hhv. q1 , q2 og q3 som vægte. Slutteligt sammenlignes dette med prisudviklingen i S&P500-indekset. Resultatet ses på Figur 4. Det bemærkes at det første
principal-aktiv overordnet set følger markedets udvikling. Det 2. principalaktiv følger nogenlunde markedets svingninger, dog mere afdæmpet. Det 3.
principal-aktiv viser klare tegn på modsatrettede bevægelser i forhold til markedet. Fx ses en kraftig nedtur omkring år 2000 efterfulgt af en kraftig optur,
hvilket er den modsatte bevægelse i forhold til markedet.
Det skal her understreges, at S&P500 løbende justeres mht. aktiver og vægte,
mens principal-aktiverne dannes som faste indekser for hele perioden. Endvidere gælder det at SP84-datasættet kun indeholder virksomheder som har
været på børsen siden 1984, og altså ikke det fulde grundlag for S&P500indekset. Der er derfor ingen grund til at forvente, at vi skulle være i stand
til at matche S&P500 fuldstændig.
Til yderligere at underbygge fortolkningen af de første par principal-aktiver
viser Figur 5 en opgørelse for SP2000 datasættet, der viser hvilke positioner,
de første 3 egenaktiver har i forskellige sektorer.
22
Figur 4: Øverst vises udviklingen i S&P indekset holdt op imod udviklingen i de første tre principal-komponenter set hver for sig. I midten
og nederst vises en vægtet sum af de tre principal-aktiver holdt op imod
S&P500. Det ses at de første tre principal-aktiver i kombination i rimelig høj grad kan forklare markedsudviklingen pånår perioden omkring
dotcom-boblen i år 2000. Vægtene er dannet vha. formlerne for efficient
rand (9).
5 Effektive beregningsmetoder
Vi vil nu kigge nærmere på hvorledes principal-komponenterne kan bruges
til at lave effektive approximationer og beregningsmetoder af historisk varians for porteføljer. Med ’effektiv’ menes en reduktion i de påkrævede antal
beregningsoperationer med en eller flere størrelsesordener, dvs. vi ønsker at
reducere beregningstiden med 90%-99%.7
For at kunne diskutere effektivitet i beregningsmetoder har vi brug for lidt
simple konventioner. Vi antager metoden implementeres i et givet programmeringssprog, der afvikles på en computer/server, hvor de basale aritmetiske
operationer (multiplikation, addition osv. ) tager lige lang tid. Vi ser bort fra
detaljer omkring repræsentation af tallene (floating points vs. fixed decimals
osv.), nøjagtighed ved udregning (altså om der accepteres afrundingsfejl) osv.
Det eneste vi er interesserede i er antallet af basale aritmetiske operationer.
Vi antager også at n >> 1, hvilket medfører at vi automatisk laver afrundinger
såsom n + 1 ≈ n og n − 1 ≈ n.
Som eksempel kan vi betragte en projektion (vektorprodukt):
ab0 = a1 b1 + a2 b2 + · · · + an bn
Det ses umiddelbart at denne udregning kræver n multiplikationer efterfulgt
af n − 1 additioner, dvs. ialt ca. 2n operationer.
6 Approximation af varians
Vi betragter igen en portefølje P med tilhørende vægte w. Vi forestiller os også, at der foreligger en tidsserie med T målinger af relative prisændringer for de
underliggende aktier, således at kovariansmatricen C kan bestemmes, og herudfra er principal-komponenterne q1 , . . . , qn og egenværdierne λ1 , . . . , λn fundet. Slutteligt er β-koefficienterne for P fundet ved at projicere vægt-vektoren
på principal-komponenterne: βi = w · q0i .
Som det fremgår af (8) er det let at finde variancen for en portefølje, når man
først har fundet principal-komponenterne, egenværdierne og beta-koefficienterne.
Dog gælder det, at udover de fortolkningsmæssige og analytiske fordele det
kan give at arbejde med principal-komponenter, vil denne fremgangsmåde ikke
i sig selv give en væsentlig beregningsmæssig besparelse af følgende årsager:
7
Indenfor datalogien ville man formulere dette som et ønske om, at reducere den beregningsmæssige kompleksitet, fx at gå fra O(n2 ) til O(n) i køretid. Denne opgave forsøger
dog at holde fokus på indsigten i PCA og de finansielle fortolkninger heraf og derfor bliver
spørgsmålet om reduktion i beregningsmæssig kompleksitet behandlet en smule forsimplet.
24
Figur 5: Øverst vises en opgørelse over positionen af de tre første
principal-aktiver indenfor forskellige sektorer opgjort for SP2000 datasættet. Positionen er opgjort ved at tælle hvor mange aktier indenfor
hver sektor principal-aktivet har hhv. en kort, en lang eller ingen
position i. Eksempelvis ses det at principal-aktiv nr. 2 er lang i 54% af
alle Consumer Services aktier. Det ses tydeligt at den 1. principal-aktiv
stort set er lig med markedet som helhed, i den forstand at det er bygget
op af en lang position af alle aktiver i markedet. Tilsvarende ses at 2.
principal aktiv er lang i cykliske aktier og kort i defensive aktier (fx
Consumer Goods og Health Care). Dette aktiv vil altså udvise særlig
stor volatilitet i hhv. højkonjunktur og lavkunjunktur.
Nederst vises et udklip fra det bagvedliggende datamateriale. Som
eksempel ses at det 2. principal-aktiv (dvs. q2 fra formlerne i afsnit
4) har en vægt på -0.065758 - dvs. en kort position - i Allstate Corp.
indenfor industrien Financials. Tilsvarende har det 3. principal-aktiv
en position på 0.00013285 i Alcoa Inc. hvilket fortolkes som ’ingen
position’. Der er anvendt den konvention at hvis en vægt ligger i
intervallet [−0.005, 0.005], da fortolkes det som hverken kort eller lang,
også kaldet ’ingen’ position.
Figur 6: Viser et tænkt men illustrativt eksempel på varianserne (dvs.
egenværdierne) af principal-aktiverne i et marked med 20 aktiver. Fx har
principal-aktiv 1 en varians på 20, principal-aktiv 2 har en varians på
10 osv. Halen af principal-aktiver har en meget lav varians omkring 1
og desuden er variansen i halen næsten konstant set imod størrelsen af
de første egenværdier.
• For at finde principal-komponenterne skal PCA-analysen først gennemføres, dvs. givet C skal man finde Q og Λ så ligning (3) er opfyldt.
Denne analyse er i sig selv krævende og baserer sig på komplekse og
beregningstunge matematiske operationer, jævnfør [5] afsnit 11. Naturligvis findes der masser af standardiserede software-pakker, som udfører
PCA-analysen, men det ændrer ikke ved at det generelt set er en beregningskrævende operation der har en løbetid, der skalerer med n2 [11].
• For at finde βi ’erne skal der laves n projektioner af vektorer med hver n
elementer i. Det vil i sig selv kræve ca. 2n2 operationer.
• Metoden forudsætter at C skal bestemmes. Ydermere gælder det at en
nøjagtig bestemmelse af alle egenværdier og principal-komponenter kræver en nøjagtig bestemmelse af C. Hvilket igen betyder at der skal anvendes en forholdsvis lang tidsserie. Hvis tidsserien indeholder T punkter
kræver udregningen af kovariansmatricen ca. T · n2 beregninger.
For at få gavn af PCA-analysen til at lave effektive beregningsmetoder af
varians m.m. er det nødvendigt at kigge nærmere på de egenværdier, som
typisk findes for finansielle tidsserier.
I praksis vil man ofte finde at der er et kraftigt fald i egenværdierne afløst af
en flad hale, et typisk billede er illustreret på Figur 6. Dette giver mulighed
for en effektiv approximation eftersom værdierne i halen kan udskiftes med en
gennemsnitsværdi.
26
Lad os anvende ideen til at regne videre på formel (8)
σ(P)2 = Var(P) = β12 λ1 + β22 λ2 + · · · + βn2 λn
2
≈ β12 λ1 + · · · + βk2 λk + λtail (βk+1
+ · · · + βn2 )
2
= β12 λ1 + · · · + βk2 λk + λtail βtail
hvor k markerer punktet lige inden egenværdi-halen starter (på figuren er k
ca. 3), λtail er den gennemsnitlige egenværdi i "halen"(på figuren ca. 1) og vi
2
2
har indført βtail
≡ βk+1
+ · · · + βn2 . Vi antager at k << n.
Hernæst skal man indse at:
2
βtail
=
n
X
βi2
i=1
−
k
X
βj2
2
= kβk −
j=1
= kw · Qk2 −
k
X
βj2
j=1
k
X
βj2 = kwk2 −
k
X
j=1
βj2
j=1
hvor det er udnytttet at Q er længde-bevarende. Slutteligt skal man se at:
(n − k) · λtail =
n
X
λi =
n
X
λi −
k
X
i=1
i=k+1
= Tr(Λ) −
k
X
λj
j=1
λj = Tr(C) −
j=1
k
X
λj
j=1
hvor det er udnyttet at Tr(Λ) = Tr(Q0 CQ) = Tr(CQQ0 ) = Tr(C).
P
Hvis vi til slut bemærker at Tr(C) = i σ 2 (Ai ) kan vi opsummere følgende
centrale resultat:
σ 2 (P) = Var(P) ≈ β12 λ1 + · · · + βk2 λk
1
+
n−k
X
σ 2 (Ai ) −
i
hvor βi ≡ w ·
k
X
j=1
!
λj
kwk2 −
k
X
!
βj2
j=1
q0i
(10)
På trods af at formlen muligvis ser lidt kompliceret ud, indeholder den væsentlige beregningsmæssig reduktioner. Bemærk især følgende:
27
• Formlen involverer kun de første k egenvektorer og egenværdier.8 De
fleste egenværdi-algoritmer fungerer iterativt på en måde så egenvektorerne findes i faldende orden.9 Det betyder også, at der opnås en væsentlig beregningsmæssig besparelse, hvis kun de første k egenværdier
og egenvektorer skal findes (husk vi antager k << n).
• Formlen involverer udover de første k egenvektorer og egenværdier kun
vektor-længden af porteføljevægtene (kwk) - som er givet på forhånd og summen af varianserne af de oprindelige aktier i markedet.
• Det kan vises teoretisk og praktisk at de første egenvektorer og egenværdier findes med højere statistisk signifikans end de efterfølgende egenvektorer og egenværdier. Ligeledes kan det vises teoretisk og praktisk at
varianser for en given tidsserie kan bestemmes med en væsentlig højere
statistisk signifikans end kovarianser. Det betyder at man ofte kan reducere antallet af målepunkter i tidsserien og stadig fastholde den ønskede
nøjagtighed i de endelige resultater.
6.1 Bestemmelse af reduktionsfaktor
Lad os nu forsøge at bestemme et kvantitativt bud på, hvilken beregningsmæssig reduktion vi opnår med formel (10). Vi antager at kovarians-matricen
C er fundet.
Først betragter vi den traditionelle fremgangsmåde: Var(P) = wCw0 . Produktet C·w0 koster n2 multiplikationer og n2 additioner, dvs. ialt 2n2 operationer. Dernæst foretages den sidste vektor-multiplikation w · (Cw0 ) som
vi allerede har set koster 2n operationer. Dvs. samlet involverer udregningen
2n2 + 2n ≈ 2n2 operationer, hvor vi antager at n er stor, dvs. n2 >> n.
Dernæst betragter vi formel (10): Vi starter med at antage at de første k egenværdier og egenvektorer er fundet. De k β-værdier involverer hver en vektorprojektion som koster 2n projektioner hver, dvs. ialt 2kn operationer. DernæstP
multipliceres med λ1 , . . . , λk , hvilket er yderligere k operationer. Hernæst
skal i σ 2 (Ai ) findes, hvilket svarer til summen af diagonalen af kovariansmatricen, dvs. ialt n operationer. Summen af de første k egenværdier koster
k operationer. Længden af w koster nP
multiplikationer og n additioner, dvs.
ialt 2n operationer. Slutteligt koster j βj2 k multiplikationer og k additioner, dvs. ialt 2k operationer. Samlet set vil udregningen af variansen ud fra
formlen således koste 2kn + k + n + k + 2n + 2k = 2kn + 4k + 3n ≈ (2k + 3)n
operationer hvor det er antaget at k << n.
Hermed får vi altså en reduktion på (2k+3)n
= k+1.5
, dvs. hvis n er ca. 10 gange
2n2
n
større end k da er reduktionen på ca. 90%. Hvis n er ca. 100 gange større end
8
9
Egenvektorerne er implicit involverede gennem β-koefficienterne eftersom βi = w · q0i .
Hvis to egenværdier ligger tæt kan de dog godt findes i omvendt rækkefølge.
28
k, da er reduktionen på ca. 99%.
6.2 Empiriske resultater
Det skal nu undersøges hvor godt approximationen (10) virker i praksis. Der
er lavet en undersøgelse for de to datasæt, SP2000 og NYSE77 og resultatet
er vist på Figur 7 og Figur 8 respektivt.
Det er vigtigt at forstå at approximationen både skal undersøges mht. forskellige valg af cutoff k og forskellige valg af porteføljer. Det giver i sig selv to
dimensioner i analysen. Dimensionen for valg af k er simpel og en-strenget.
Dimensionen for valg af portefølje er derimod kompleks. Forskellige typer af
porteføljer virker nemlig meget forskelligt på formlen. Eksempelvis finder man
at diversificerede porteføljer hvor mange aktier indgår, egner sig langt bedre
til at blive approximeret end porteføljer som består af et lille antal aktier, fx.
4. Lad os vende tilbage til dette lidt senere i dette afsnit. Først gennemgås
fremgangsmåden i analysen:
• Der udvælges M porteføljer, hvor M er valgt til at være lig n, men det
kunne lige så godt have været et andet antal såsom 500 eller 1000. Se
nedenfor hvordan disse porteføljer er udvalgt.
• For hver af de M porteføljer og for hver cutoff-værdi k = 1, 2, . . . , n
udregnes den approximerede standardafvigelse via formel (10).
• Samtidig udregnes for hver portefølje den historiske - dvs. den sande standardafvigelse.
• Herefter findes for hvert k og hver portefølje den procentvise fejl mellem
den sande standardafvigelse og den approximerede, som herefter plottes
med k ud af x-aksen og den procentvise fejl på y-aksen. Disse plot ses
på Figur 7 for SP2000 og Figur 8 for NYSE77-datasættet. Bemærk at
for hver k-værdi haves et målepunkt for hver af de M porteføljer hvilket
giver anledning til et tykt ’bånd’ af grafer. Tætheden i grafen viser hvor
mange porteføljer, der havde en given approximationsfejl.
NB: Det skal bemærkes at vi forsøger at gennemføre en analyse med nogle
statistiske størrelser, som nemt kan forvirre, fordi vi så at sige forsøger at finde
en standardafvigelse af fejlen på en approximation af en standardafvigelse.
Derfor skal man holde tungen lige i munden.
6.2.1 Historisk std.afv. for en portefølje
Lad os kort redegøre for, hvorledes den historiske std.afv. findes for en konkret
portefølje. Husk på at en portefølje er repræsenteret ved en vektor w med
29
Figur 7: Viser den relative fejl på standardafvigelsen udregnet efter formel (10) som funktion af cutoff-faktor for datasættet SP2000.
Figur 8: Viser den relative fejl på standardafvigelsen udregnet efter formel (10) som funktion af cutoff-faktor for datasættet NYSE77.
vægte, der summer til 1, hvor wi angiver hvor stor en procentdel af porteføljens
værdi, som udgøres af den i’te aktie i markedet.
Husk også på at de grundlæggende markedsdata består af en tidsserie som vi
kan betegne med G (som står for ’gain’). Tidsserien angiver for hver dag det
relative afkast for hver enkelt aktie, og vi tænker på G som en matrix hvor
Gij angiver afkast på tidspunkt i for aktie j.
Nu kan vi projicere afkastene for aktierne i markedet ned på porteføljevægtene
ved at udregne matrixproduktet g0 ≡ G·w0 . Hermed finder man en ny tidsserie
i form af en søjlevektor, som repræsenterer de daglige relative afkast man ville
have haft, hvis man havde haft den pågældende portefølje. Slutteligt finder
man den historiske standardafvigelse på sædvanligvis fra denne tidsserie. Som
√
matrixprodukt ser det således ud: σ = g · g0.
6.2.2 Valg af porteføljer
Når porteføljerne skal vælges til analyse af (10) er det værd at have en overordnet ide om, hvilke typer af porteføljer der er interessante at betragte. Det
viser sig at den vigtigste ’parameter’ i forhold til om porteføljen egner sig til
approximation eller ej, er hvorvidt porteføljen er diversificeret eller ej.
Forklaringen er at diversificerede porteføljer naturligt har en stor vægt på de
første principal-aktiver og samtidig vægte spredt godt ud i halen, mens ikkediversificerede porteføljer kun har lille vægt på de første principal-aktiver og
til gengæld har stor vægt på enkelte principal-aktiver i halen. Det betyder at
approximationen bliver dårlig, jævnfør Figur 6. Hvis det ikke umiddelbart er
klart ud fra denne forklaring, kan vi gå lidt videre med nogle observationer
omkring Figur 6.
Læg først mærke til at de første egenværdier (og vi minder om at egenværdierne eksakt er lig variansen af principal-aktiverne) er flere størrelsesordener
større end dem i halen. Det betyder hvis blot man har en stor vægt på de første principal-akter, da vil bidraget fra halen være af mindre betydning. Læg
hernæst mærke til, at halen har en svag hældning. Faktisk er det oftest sådan
at halen falder nogenlunde linært ned til 0, eller i hvert fald meget tæt på
0. Husk nu på at kernen i approximationen går ud på at udskifte de eksakte
værdier i halen med en gennemsnitlig værdi for halen. Det betyder at hvis en
given portefølje har en stor vægt midt i halen, da vil approximationen gå godt
(fordi gennemsnitsværdien af halen er meget tæt på den midterste værdi af
halen). Men hvis en given portefølje har en stor vægt i en af siderne i halen, da
vil approximationen gå meget galt. Til sidst bemærkes at hvis en portefølje
har lige meget vægt i begge sider af halen, da vil approximationen igen gå
godt.
For at lave en undersøgelse af approximationens gyldighed for porteføljer med
forskellige grader af diversifikation, genererer vi et antal porteføljer med én
32
aktie i hver, et antal porteføljer med 4 aktier i hver, osv. med 8, 16, 32 og 64
aktier i hver portefølje. Til sidst laver vi et antal porteføljer hvor alle aktier
indgår med en tilfældig valgt vægt, hvorefter der normeres så summen bliver
1.
Tabel 2: Cutoff-værdier for approximation med fejl under 1% vist for
forskellige grader af aktie-diversifikation angivet med 1, 4, 8, ..., 64,
’Alle’ samt for forskellige datasæt: NYSE77, SP2000 og SP84. Øverste
tabel viser de absolutte cutoff-værdier, nederste tabel viser cutoff-værdien
relativt til det totale antal principal-aktiver/egenværdier.
Cutoff
NYSE77
SP2000
SP84
Cutoff %
NYSE77
SP2000
SP84
1
383
109
53
1
22%
29%
54%
4
307
63
38
4
18%
17%
39%
8
174
13
2
8
10%
3.4%
2.0%
16
32
64
Alle
29
1
1
1
1
1
1
1
1
1
1
1
16
32
64
Alle
1.7% 0-1% 0-1% 0-1%
0-1% 0-1% 0-1% 0-1%
0-1% 0-1% 0-1% 0-1%
6.2.3 Resultat af empiriske undersøgelser
Resultaterne af de empiriske undersøgelser af approximationsformlen (10) for
standard-afvigelsen er opsummeret i Figur 7, Figur 8 og Tabel 2.
Graferne på de to figurer viser hvilken procentvis fejl der findes for de valgte
porteføljer som funktion af cutoff-factor k. I den øverste graf på Figur 7 kan
man eksempelvis se at approximation af std.afv. for en portefølje med kun
én aktie vil give fejl på mellem 0 og 40%, hvis k = 1 (dvs. kun den første
egenværdi medregnes). Hvis k = 200 svinger fejlen mellem 0 og godt 10%.
Hvis alle egenværdier medregnes (dvs. k er lig n som er lig 381) er fejlen
eksakt 0, som forventet.
Hvis vi går til næste graf for porteføljer med 4 aktier i hver, ser man det
samme mønster, dog er approximationen væsentlig bedre allerede når k er
tæt på 1. Hvis man kan leve med en præcision på 5% kan man faktisk nøjes
med at sætte k = 5, hvilket i sig selv vil give en reduktion i beregningstiden
på 1 − 5/381 ≈ 98.7%.
For porteføljer med 16 aktier og opefter kommer fejlen under 1% allerede med
k = 1.
Som man kan se er der meget stor forskel på hvor godt den historiske varians for forskellige porteføljer egner sig til at blive approximeret, hvilket giver
anledning til et bredt bånd af grafer. Til gengæld ser man at diversificerede porteføljer generelt set danner et smallere bånd end ikke-diversificerede
33
porteføljer og at diversificerede porteføljer egner sig væsentlig bedre til approximationen.
På Tabel 2 ses øverst hvilken cutoff-værdi k man skal bruge for at få en
approximationsfejl på under 1% for porteføljer med forskellig grad af diversifikation. Fx. kan man se at for datasættet NYSE77 og en portefølje med 16
aktier, skal man bruge de første 29 principal-aktiver, dvs. k = 29, for at få en
approximationsfejl på under 1%. Tilsvarende kan man aflæse at med en fuldt
diversificeret portefølje kan man nøjes med det første principal-aktiv, dvs.
k = 1. Den nederste tabel viser det samme, blot er tallene omregnet til en reduktionsfaktor, dvs. man finder cutoff-faktoren divideret med det totale antal
principal-aktiver: k/n. Dette er samtidig et udtryk for den mulige reduktion
man kan opnå i beregningstiden, jævnfør afsnit 6.1. Eksempelvis kan man for
porteføljer med mindst 16 forskellige aktier i SP2000 (som er et marked med
384 aktier) få en reduktion på ca. 100% − 3.4% = 96.6% af beregningstiden
og stadig have en præcision på 1%.
Ved generelt at kigge på graferne for diversificerede porteføljer kan man se, at
approximationen virker med meget stor nøjagtighed.
7 Definition af value-at-risk og expected shortfall
Vi har nu vist, hvordan man kan finde en god approximation af den historiske
varians og standard-afvigelse for en givet portefølje. I praksis er man dog ofte
mere interesseret i andre risikomål end standard-afvigelsen. Typisk ønsker
man nemlig risikomål som kan bruges til at bestemme hvor meget risikovillig
kapital en given position i markedet kræver, hvilket fører til en betragtning
af risikomål som giver et bud på hvor store tab, man kan blive udsat for.
To sådanne risikomål er value-at-risk (VaR) og expected shortfall (ESF) som
behandles i dette afsnit.
VaR er løst sagt et mål for hvor stort et værditab vi højst kan udsættes for
indenfor en given tidshorisont og med en given sandsynlighed. For at definere
VaR mere formelt tager vi udgangspunkt i en given portefølje P, et givet
tidsinterval t - som i rammerne af denne analyse altid er 1 dag - og en given
sandsynlighed α ∈ [0, 1], typisk 95% eller 99%. Det er praktisk at indføre α
ˆ=
1 − α, dvs. α
ˆ vil typisk være 5% eller 1%. Nu minder vi om at α-fraktilen for
P angiver det punkt på x-aksen i et diagram over sandsynlighedsfordelingen
for P, hvor α-% af sandsynlighedsmassen ligger til venstre for punktet og
dermed ligger 1 − α-% til højre for punktet.
Hermed kan vi definere value-at-risk sådan her:
VaRα (P) = −1 · (ˆ
α-fraktilen af P)
(11)
Det negative fortegn er indført, fordi man ønsker at value-at-risk skal være et
34
Figur 9: Illustration af sandsynlighedsfordeling og fraktiler. Figuren viser et histogram over daglige relative afkast for Alcoa Inc. (med ticker
symbolet AA) målt over perioden 1984-2014, hvilket svarer til en tidsserie med ca. 7500 målepunkter. Nederste billede viser halerne i større
detalje. Middelværdien kan findes til at være 0.0004 og spredningen er
0.023. Ved hjælp af lidt databehandling kan man finde, at 5%-fraktilen
ligger omkring -0.033, dvs. med andre ord er VaR95% = 0.033. Dernæst
kan man finde at 1%-fraktilen ligger omkring -0.06, dvs. VaR99% = 0.06.
NB: Hvis det daglige afkast havde været normalfordelt, skulle 99.7% af
punkter ligge mellem -0.07 og +0.07, hvilket passer nogenlunde grafisk
set. Til gengæld burde et daglig relativt afkast på -0.24 (som er det yderste venstre punkt) kun ske en ud af ca. 8 · 1026 dage, hvilket er en klar
indikation af, at det daglige afkast ikke er normalfordelt.
Figur 10: Grafisk illustration af VaR og ESF. Her vises som eksempel
90%-VaR og som det fremgår ligger 90% af sandsynlighedsmassen til
højre for punktet mens 10% ligger til venstre. ESF kan man forstå som
middelværdi af den del af halen som ligger til venstre for VaR, dvs.
middelværdien af den gule hale.
positivt mål, dvs. det angiver det potentielle tab. Begreberne omkring fraktiler
og definitionen af VaR er illustreret på Figur 9 og Figur 10.
NB: Bemærk at value-at-risk traditionelt defineres som et absolut mål for det
potentielle tab og value-at-risk måles normalt i en konkret valuta. Fx kan
man have en 95% VaR på 200.000 kr for en portefølje med værdi 1mio.kr.
Som allerede nævnt tænker vi på P som en stokastisk variabel, der angiver
den daglige relative tilvækst. Hvis fx P en dag er lig -0,08, da svarer det til
at porteføljen har haft et relativt afkast på -8%, eller omvendt sagt, et tab på
+8%. Da vi benytter den konvention, at alle porteføljer koster 1 kr (og bliver
rebalanceret hver dag) kan vi lige så godt tænke på det som om, at vi har lidt
et tab på 0,08kr. Med andre ord kan vi tænke på vores stokastiske variabel P
som den daglige gevinst/tab i kroner, og dermed bliver ovenstående definition
af value-at-risk konsistent med den normale definition.
Expected shortfall (forkortet ESF) er løst sagt et mål for, hvor meget man
skal forvente at tabe, hvis det først går galt. Mere præcist defineres ESF i
lighed med VaR for en given sandsynlighed α% og en given tidshorisont, som
det forventede tab i de α-% værste scenarier der kan hænde indenfor den givne
tidshorisont. ESF er også illustreret på Figur 10.
Formelt kan vi definere ESF som følger:
ESFα = −E(P|P < −VaRα )
Igen er der indført et negativt fortegn, for at ESF bliver en positiv størrelse. I
ord vil man læse formlen som det forventede tab givet et tab større end VaRα .
36
7.1 Historisk VaR og ESF
Når man arbejder med tidsserier for daglige relative prisændringer er det let at
finde den historiske VaR og historiske ESF. Historisk VaR for en givet fraktil
α findes ved at sortere tidsserien i stigende orden efter daglig relativ tilvækst
(så de største tab kommer først) og dernæst finde det K’te element i rækken,
hvor K = α
ˆ T og T er det totale antal elementer i tidsserien. Vi minder om at
α
ˆ ≡ 1 − α.
Historisk ESF findes i forlængelse heraf, ved at tage middelværdien af de
elementer, som kommer før det K’te element.
Det skal her bemærkes at historisk VaR og ESF også kan bestemmes på andre
måder fra en given tidsserie. Fx kan der benyttes en form for interpolation,
hvis man forsøger at bestemme nogle fraktiler, hvor der er langt imellem punkterne. Dette behov opstår typisk hvis man har et lille antal datapunkter i sin
tidsserie, eller hvis man forsøger at bestemme de mest ekstreme fraktiler såsom 99.9999%-fraktilen. Til denne opgave vil vi nøjes med den simple metode
beskrevet herover.
7.2 Empiriske resultater
Vi kan starte med at undersøge VaR for enkelte aktier, hvor vi sammenligner
historisk VaR med 1.6449×σ som teoretisk skulle give det samme, hvis afkastet
var normalfordelt. Dette er vist for SP1984 på Figur 11. Det ses at der som
forventet er en klar sammenhæng mellem σ og VaR i den forstand at de
bevæger sig op og ned sammen (dvs. de er stærkt korrelerede), men det ses også
at de ikke ligger oven i hinanden (dvs. afkastet er ikke normalfordelt) og der
er heller ikke perfekt korrelation, dvs. der er lidt forskel på hvor tæt de ligger
på hinanden. Dette skal fortolkes således, at der er forskelle i haletykkelserne
på forskellige aktier. Til at illustrere dette kigger vi nærmere på to specifikke
aktier, WMB (Williams Companies, Inc.) og PNR (Pentair plc). Tabel 3 viser
VaR-værdier for forskellige α-værdier. Det ses tydeligt at WMB er tykkere i
halen end PNR idet alle fraktiler udover 99%-fraktilen ligger længere ude.
Derudover kunne det være interessant at sammenligne VaR og ESF for aktier
versus principal-aktiver. Dette er vist på Figur 12. Som det kan ses giver
principal-aktiverne en form for koncentration af volatiliteten således at de
første principal-aktiver har en høj VaR svarende til en høj volatilitet, mens
halen har en lav VaR. Det samme gør sig gældende for ESF.
37
Alpha
95%
98%
99%
99.7%
99.9%
99.99%
VaRα /σ
PNR Norm.ford.
1.5
1.6
2.1
2.1
2.6
2.3
3.8
2.8
5.1
3.1
11.6
3.7
WMB
1.1
1.8
2.5
4.1
5.8
19.8
Tabel 3: Viser value-at-risk målt i enheder af σ for to virksomheder PNR
(Pentair Plc) og WMB (Williams Companies Inc.) med forskellig haletykkelse sammenlignet med en normalfordeling. Eksempelvis er VaR-95%
for Pentair Plc ca. lig med 1.5 standardafvigelser.
Figur 11: Illustration af VaR fundet henholdsvis historisk og ud fra formlen 1.6449 × σ.
Figur 12: Illustration af VaR (øverst) og ESF (nederst) fundet for alle
aktierne i SP84 markedet sammenlignet med VaR fundet for principalaktiverne. Den ’glatte’ grønne kurve angiver VaR hhv. ESF for principalaktiverne. Den blå hakkede kurve angiver VaR hhv. ESF for aktierne i
markedet.
Figur 13: En normalfordeling er entydigt kendetegnet ved dens middelværdi og standardafvigelse.
8 Approximation af value-at-risk
Vi ønsker at finde en metode til udregning af VaR og ESF der ligger i forlængelse af den allerede udledte formel (10) for approximation af varians og
standardafvigelse. Samtidig ønsker vi at forstå de begrænsninger og antagelser, der skal gøres for at approximationen er god.
Først gøres den betragtning, at visse familier af stokastiske variable har et fast
forhold mellem standard-afvigelsen og en given fraktil. For disse stokastiske
variable vil der også være et fast forhold mellem standard-afvigelsen og VaR-α,
jævnfør definitionen af VaR (11). Som eksempel har man den velkendte 68-9599.7 regel for normalfordelingen, der siger, at ca. 68% af sandsynlighedsmassen
ligger inden for 1 σ, 95% ligger inden for 2σ og 99.7% ligger inden for 3σ.
Denne ide er illustreret på Figur 13. Således har man at VaR-95% for en
normalfordelt stokastisk variabel med varians σ altid vil være 1.6449 × σ.
Man kan nu opstille den hypotese, at de porteføljer P vi kan konstruere i markedet, vil have denne egenskab, altså at der givet α er et fast forhold mellem
σ(P) og VaRα (P), uanset hvilken portefølje P vi vælger. Teoretisk set er der
gode eksempler på multivariate fordelinger, som har denne egenskab, specielt
skal nævnes de elliptiske fordelinger, jf. [2] afsnit 3.3. En simpel måde at efterprøve denne tese empirisk er, at udregne den historiske value-at-risk samt den
historiske standard-afvigelse for en række forskellige porteføljer eller aktier
i markedet, og dernæst finde forholdet mellem VaR og standard-afvigelsen.
Dette er undersøgt for SP84 og resultatet kan ses på Figur 14 og Tabel 4.
Undersøgelsen viser (som forventet) at hypotesen ikke holder eksakt. Til gengæld kan vi se, at forholdet ligger i et fast bånd, der i værste fald svinger
med ±10% for aktierne og i bedste fald ca. ±4% for halen af egenvektorerne,
jævnfør Tabel 4.
Vi kan nu lave følgende approximation, som analytisk set er grov og klodset
men som til gengæld er nem at forstå. Betragt en portefølje P og en fast værdi
40
Figur 14: Viser en undersøgelse af forholdet mellem value-at-risk og
std.afv. for SP84-aktierne øverst og SP84-principalaktiverne (dvs. egenvektorerne) nederst. Først observeres at forholdet generelt er mindre for
aktierne (svinger mellem 1.2 og 1.5) end for egenvektorerne (svinger
mellem 1.4 og 1.6). Det kan fortolkes således at egenvektorerne repræsenterer de kombinationer af aktiver i markedet, som giver de tungeste
haler! Det er ønskværdigt set ud fra et risikostyringsperspektiv, idet vi
jo netop forsøger at få styr på den risiko, som gemmer sig i de tunge haler. Dernæst ses det, at grafen for egenvektorerne har den særlige
egenskab, at halen stabiliserer sig en smule ved at svinge mellem 1.5 og
1.6. Tabel 4 viser middelværdi og std.afv. for svingningerne for hhv. aktierne, egenvektorerne og halen af egenvektorer (hvor de første 10 hhv.
20 egenvektorer skæres fra).
Tabel 4: Middelværdi og std.afv. i forholdet value-at-risk / std.afv.
Middelværdi
Aktier
1.42
Egenvektorer
1.57
Egenvektor > 10
1.58
Egenvektor > 20
1.59
std.afv.
0.071
0.055
0.035
0.03
95% konfidensinterval
±10%
±7%
±4.4%
±3.8%
af α og sæt
γ ≡ VaR(P)/σ(P)
Vi har droppet α-fodtegnet for at holde notationen simpel, men det er underforstået at vi mener VaRα () når vi skriver VaR(). Vi tænker på γ som
det sande men ukendte forhold mellem value-at-risk og std.afv. for P. Definér
ligeledes:
γi ≡ VaR(Qi )/σ(Qi )
hvor vi husker at Qi er den i’te principal komponent. Vi tænker på γi som forholdet mellem VaR og std.afv. for det i’te principalaktiv. Med udgangspunkt
i undersøgelsen fra Figur 14 og Tabel 4 vil vi nu antage at γ ≈ γi for alle i.
Ved at anvende egenvektor-dekompositionen fra (8) ser vi følgende:
(VaR(P))2 = γ 2 · σ 2 (P)
= γ 2 β12 λ1 + · · · + βn2 λn
= β12 γ 2 λ1 + · · · + βn2 γ 2 λn
≈ β12 γ12 λ1 + · · · + βn2 γn2 λn
= β12 VaR(Q1 )2 + · · · + βn2 VaR(Qn )2
hvor det i sidste linie er udnyttet at γi2 λi = γi2 σ 2 (Qi ) = (γi σ(Qi ))2 = VaR2 (Qi ).
Hvis vi indsætter α-fodtegnet igen og tager kvadratroden på begge sider finder
vi den endelige formel:
q
VaRα (P) ≈ β12 VaR2α (Q1 ) + · · · + βn2 VaR2α (Qn ) hvor βi ≡ w · q0i
(12)
Selvom argumentet i udledningen ikke er specielt stærkt, skal vi senere se at
denne approximation empirisk set er overraskende god.
Som en sidste bemærkning vil vi sammenligne denne formel med en tilsvarende
formel for standard-afvigelsen. Formel (8) kan omskrives til:
q
σ(P) = β12 σ 2 (Q1 ) + · · · + βn2 σ 2 (Qn )
(13)
42
Det understreges at σ er et moment-baseret risikomål, og som sådan er det
fundamentalt anderledes end VaR, ESF og andre fraktil-baserede risikomål.
Ligeledes skal det understreges at (13) er et analytisk eksakt resultat og ikke en
approximation som formel (12). På trods af forskellene i σ og VaR og måden,
hvorpå de to formler er udledt, så giver det intuitivt mening at riskomål som i
en eller anden forstand "skalerer"med std.afvigelsen, kan dekomponeres vha.
principal-komponenterne på samme måde som std.afvigelsen.
Ovenstående approximation har stadig den ulempe, at den involverer alle egenvektorer, hvilket både er beregningstungt og udsat for stor statistisk støj. Derfor vil vi nu forbedre approximationen yderligere ved at indsætte en middelværdi for value-at-risk i halen af egenvektorer og samtidig udskifte de mange
βi -værdier for halen med en samlet projektion ned i underrummet udspændt
af hale-egenvektorerne - præcis som det blev gjort i forbindelse med udledning
af formel (10). Hermed ender vi med dette centrale resultat, som gør det muligt at approximere value-at-risk med de k første egenvektorer og egenværdier
til kovarians-matricen:
(VaR(P))2 ≈ β12 VaR2 (Q1 ) + · · · + βk2 VaR2 (Qk )
!
k
X
2
2
2
+ VaRtail · kwk −
βj
(14)
j=1
hvor βi ≡ w ·
q0i
8.1 Bestemmelse af VaR i halen
Approximationen (14) involverer kun de første k egenvektorer bortset fra den
gennemsnitlige value-at-risk for halen af egenvektorer. I tilfældet med approximationen (10) kunne man anvende et analytisk eksakt resultat om at
summen af egenværdierne er lig summen af varianserne for de oprindelige ak2
tier (matematisk Tr(C) = Tr(Λ)) til at omskrive σtail
yderligere uden tab af
nøjagtighed. Denne mulighed har vi ikke nu. Man kan istedet forestille sig en
af følgende metoder anvendt til at bestemme VaR2tail :
1. Man kan finde alle principal-aktiverne i halen og dernæst finde value-atrisk for hver af disse. Slutteligt kan man finde gennemsnittet af VaR2i .
2. Som det kan ses på Figur 12 er halen nogenlunde jævnt faldende. Dvs.
man kunne nøjes med at bestemme value-at-risk for et par udvalgte
principal-aktiver i halen, og dernæst lave en simpel lineær interpolation
for at finde den gennemsnitlige værdi af VaR2i for i > k.
43
Uanset hvilken model der vælges til at estimere den gennemsnitlige værdi af
VaR2i i halen skal det bemærkes, at det ligesom for approximation af standardafvigelsen/variansen jf. (10) gælder, at de første principal-aktiver vil veje
tungest for vel-diversificerede porteføljer, både fordi VaR er størst for de første
principal-aktiver jf. Figur 12 og fordi vel-diversificerede porteføljer vil have de
største βi -værdier for de første i’er. Derfor vil selv større unøjagtigheder ved
bestemmelse af Var2tail blive mindre betydende når vi bruger approximationen
(12).
Til de empiriske undersøgelser i denne opgave vil vil blot anvende metode 1.
Dette kan lade sig gøre fordi vi betragter et relativt lille marked og dermed
kan vi nemt finde alle principal-aktiverne i halen.
8.2 Empiriske resultater
For at efterprøve formel (14) er der i genereret 5000 diversificerede porteføljer.
For hver af disse er VaR udregnet historisk og ved hjælp af formel (14) for forskellige cutoff-værdier k. Dernæst er den procentvise fejl mellem den historiske
og den approximerede værdi fundet og slutteligt er der lavet et histogram, der
viser fordelingen af procentvise fejl. Derudover er den gennemsnitlige fejl og
spredningen på fejlen fundet.
Undersøgelsen er først lavet uden cutoff (dvs. alle principal-aktiverne bruges)
på tværs af datasættene SP84 og SP2000. Resultatet ses på Figur 15.
Det ses at approximationen er brugbar, hvis man kan leve med unøjagtigheder på VaR på op til ca. 4%. Bemærk at denne fejlmargin opstår pga. den
approximation der blev lavet ved udledning af formel (12) og ikke pga. haleapproximationen i formel (14), da den ikke er anvendt her.
For at se om hale-approximationen giver en yderligere forringelse af nøjagtigheden, laves en undersøgelse af SP2000, hvor der prøves med forskellige
cutoff-værdier, denne gang for VaR-99% istedet for VaR-95%. Resultatet ses
på Figur 16. Det bemærkes at alle 4 histogrammer har en forholdsvis høj grad
af overensstemmelse hvad angår spredning og form. Det skal fortolkes således,
at man ved at anvende et meget lille principal-aktiver i approximationen kan
bestemme VaR med nogenlunde samme præcision, som hvis man anvender
alle principal-aktiver.
Approximationen af VaR er måske ikke fuldt tilfredsstillende, idet en fejlmargin på 4-6% ikke nødvendigvis er acceptabel. Da udregning af risikomål typisk
handler om kapitalkrav, er det ikke ligegyldigt, om der skal bruges 4-6% mere
eller mindre. Vi skal nu se, at ESF som risikomål har nogle pænere egenskaber
i denne sammenhæng, da vores approximation virker væsentlig bedre for ESF
end for VaR.
44
Figur 15: Viser et histogram over den procentvise fejl mellem den historiske og den approximerede værdi af VaR-95% jf. formel (14) for hhv.
SP84 (øverst) og SP2000 (nederst) for 5000 tilfældigt genererede porteføljer.
Figur 16: Viser den procentvise fejl mellem den historiske og den approximerede værdi af VaR-99% for SP2000 med anvendelse af cutoffværdier på hhv. 1, 4 og 8 (de tre øverste) mens det nederste histogram
er lavet uden cutoff.
9 Approximation af expected shortfall
Vi vil nu gå skridtet videre og lave den samme approximation for expected
shortfall. Igen gør vi den observation at der for visse familier af stokastiske variable vil være et fast forhold mellem standardafvigelsen og expected shortfall
for en given fraktil. Dette vil eksempelvis gælde for de multivariate elliptiske
fordelinger.
Ved at bruge de samme argumenter som blev anvendt til udledning af (14),
finder man en approximation for ESF som følger:
ESFα (P) ≈
q
β12 ESF2α (Q1 ) + · · · + βn2 ESF2α (Qn ) hvor βi ≡ w · q0i
(15)
Lige som før kan vi også anvende argumentet med at bruge en gennemsnitsværdi for ESF2α (Qi ) i halen for i > k, dvs. vi finder formlen:
(ESF(P))2 ≈ β12 ESF2 (Q1 ) + · · · + βk2 ESF2 (Qk )
!
k
X
2
2
2
+ ESFtail · kwk −
βj
(16)
j=1
hvor βi ≡ w ·
q0i
Ligesom for approximationen af VaR skal vi overveje, hvordan vi kan finde
ESF2tail , dvs. den gennemsnitlige expected shortfall for principal-aktiverne i
halen. Og ligesom for VaR kan man forestille sig forskellige heuristiske metoder
til at finde dette, jævnfør diskussionen i afsnit 8.1. Til denne opgave bruger vi
samme fremgangsmåde for ESF som vi gjorde for VaR, dvs. vi finder ganske
simpelt expected shortfall for alle principal-aktiverne med index i>k og herfra
finder vi den gennemsnitlige værdi af expected shortfall.
9.1 Empiriske resultater
Efterprøvningen af approximationerne for ESF foregår efter samme metode
som for VaR og resultaterne ses på Figur 17, Figur 18 og Figur 19.
Overordnet set bemærker vi at approximationerne generelt er bedre for ESF
end for VaR. Fx. ses på Figur 17 at approximationen har en fejlmargin på
under 1% for både SP84 og SP2000.
47
Figur 17: Viser den procentvise fejl mellem den historiske og den approximerede værdi af ESF-95% for hhv. SP1984 (øverst) og SP2000 (nederst).
Figur 18: Viser den procentvise fejl mellem den historiske og den approximerede værdi af ESF for SP2000 for forskellige værdier af α lig
90%, 95%, 99% og 99.9% (øverst til nederst). Det ses at approximationen
er ret god for alle α-værdier op til og med 99%, hvorefter approximationen bliver forringet.
Figur 19: Viser den procentvise fejl mellem den historiske og den approximerede værdi af ESF for forskellige cutoff-værdier k = 1, 4, 8 for
SP2000. Den nederste figur er lavet uden cutoff.
Fra Figur 18 kan vi observere at fejlmarginen vokser når fraktilen α vokser.
Dvs. vores approximation bliver dårligere og dårligere jo længere vi kommer
ud i halerne. For at forstå dette fænomen, skal vi tænke lidt dybere over,
hvad der gør principal-aktiverne så velegnede til bestemmelse af variansen,
faktisk så velegnede at formel (8) gælder eksakt. Forklaringen er, at principalaktiverne udgør de aktiver i vores marked, som maksimerer variansen10 og
som samtidig er ukorrelerede. Hvis vi samtidig husker at variansen og korrelation pr. definition er 2. ordens momenter, kan vi omformulere og sige, at
principalaktiverne er den optimale basis til at beskrive og beregne 2. ordens
momenter i markedet. Samtidig skal man indse at det ikke nødvendigvis er
optimalt til at beskrive og beregne højere ordens momenter.
Da VaRα og ESFα er fraktil-bestemte mål og α desuden kan antage alle værdier
mellem 0 og 1, afhænger de ikke kun af ét moment, men derimod af alle
momenter.11 Den grundlæggende hypotese for denne opgave kan nu udtrykkes
således:
Principalaktiverne udgør en optimal basis til at beskrive og beregne varians og standard afvigelse. Samtidig udgør de et godt bud
på en tæt-på-optimal basis til at beskrive og beregne VaR og ESF.
10
Og her er det underforstået at vi mener, at den første egenvektor maksimerer variansen
globalt, den anden maksimerer variansen i underrummet vinkeltret på den første, den tredje
maksimerer variansen i underrummet vinkelret på de to første osv. osv.
11
Vi minder om at en sandsynlighedsfordeling (for ’pæne’ fordelinger) entydigt bestemmes
af sine momenter og omvendt, jævnfør Teorem 5.2 i [4] som vedrører moment-genererende
funktioner.
51
10 Bestemmelse af cut-off faktor vha. af random matrix
teori (RMT)
Random matrix teori har sit udspring i teoretisk fysisk hvor grundstenene
blev lagt af Eugene Wigner, som i øvrigt vandt Nobelprisen i fysik i 1963.
Random matrix-teori har fundet anvendelse indenfor en lang række forskellige
områder indenfor fysisk, matematik og statistik og i de senere år også indenfor
finansiering, jævnfør [6] og [7].
I random matrix-teori betragter man ganske simpelt en matrix, hvor de enkelte
elementer er stokastiske variable eller - om man vil - tilfældige tal udtrukket
efter en given sandsynlighedsfordeling. Herefter undersøger man forskellige
statistiske egenskaber ved matricen, som udgangspunkt relateret til fordelingen af egenværdier og egenvektorer.
Vi vil i det følgende betragte et marked med n aktiver, som antages at være
ukorrelerede og med identiske sandsynlighedsfordelinger. Bemærk at vi ikke
antager at de er normalfordelte, blot at de er identisk fordelte. For en nemheds
skyld antager vi også at de er normaliseret så standardafvigelsen er 1. Det
antages nu at vi har en tidsserie i dette marked der viser afkast (pr. minut,
pr. time, pr. dag eller andet) af længde T , som angiver antal målepunkter.
Følgende må gælde helt oplagt:
• Den ægte (dvs. teoretiske) korrelationsmatrix for markedet må være lig
identitetsmatricen:


1 0 ... 0
0 1
0
.

.

C teoretisk = I n = ..
. . ... 
0
0
...
1
• Identitetsmatricen I n ses at have n egenværdier alle med værdien 1 mens
egenvektorerne er givet ved enhedsvektorerne.
• Når T bliver meget stor for fastholdt n må den empiriske korrelationsmatrix nærme sig den teoretiske grænse, dvs. fordelingen af egenværdier
må konvergere mod en spike, som har meget stor tæthed omkring 1.
• Jo mindre T bliver (dvs. jo færre målepunkter der haves i tidsserien) jo
større usikkerhed må der være i bestemmelsen af den empiriske korrelationsmatrix, dvs. naivt set må man forvente at fordelingen af egenværdier spreder sig i en eller anden form omkring 1.
Et centralt resultat indenfor RMT er, at der kan findes en overraskende simpel
analytisk formel der beskriver fordelingen af egenværdier i grænsen hvor både
52
Figur 20: Viser den teoretiske grænse for distributionen af egenværdier for korrelationsmatricen for en random matrix for Q ≡ T /N =
1.5, 10, 100 respektivt, hvor Q = 1.5-værdier svarer til den brede fordeling med en peak i venstre side mens Q = 100 svarer til den højeste
og mest centrerede halv-cirkel. Det man skal lære af figuren er, at jo
kortere tidsserie man har relativt til markedets størrelse (dvs. jo mindre
Q-værdi), jo mere bredspektret bliver støjen. Omvendt, jo længere tidsserie man vælger, jo mere konvergerer formen mod en spike omkring 1,
der i den teoretiske grænse samler hele sandsynlighedsmassen præcis i
punktet 1.
n og T går mod uendelig på en sådan måde at forholdet Q ≡ T /n fastholdes.
Hvis vi lader ρ(λ) betegne tæthedsfunktionen for fordelingen af egenværdier
ser formlen således ud:
p
Q (λmax − λ)(λmin − λ)
ρ(λ) =
λ
2π
p
max
hvor λmin = 1 + 1/Q ± 2 1/Q og Q = T /n
(17)
Formlen er oprindelig udledt af Vladimir Marcenko og Leonid Pastur i 1967, og
fordelingsfunktionen kaldes derfor for en Marcenko-Pastur fordeling. Fordelingen kan opfattes som et støjbånd, dvs. egenværdier som falder indenfor båndet
har ingen statistisk signifikans. Marcenko-Pastur fordelingen er illustreret på
Figur 20 for tre forskellige Q-værdier.
For at efterprøve disse resultater, kan vi simulere en tidsserie af længde T
fra et fiktivt marked bestående af n ukorrelerede aktiver med daglige afkast
der følger fx en normalfordeling (men vi kunne have valgt en vilkårlig anden
fordeling med middelværdi 0 og spredning 1). Herudfra findes korrelationsmatricens egenværdier og der kan tegnes et histogram. Dette sammenlignes med
Marcenko-Pastur fordelingen. Resultatet ses på Figur 21 hvor sammenfaldet
mellem empiri og teori er tydeligt.
53
Figur 21: Empirisk egenværdi-fordeling af en korrelationsmatrix for et
fiktivt marked bestående af ukorrelerede og identisk fordelte aktiver hvor
markedet størrelse n og tidsseriens længde T er valgt til at matche hhv.
SP2000 (øverst) og SP84 (nederst). Til sammenligning er også indtegnet
den teoretiske Marcenko-Pastur fordeling (tynd kurve). Det ses at den
empiriske fordeling tilnærmes fint af den teoretiske. Det skal understreges at i grænsen hvor både n og T går mod uendelig i et fast forhold,
da vil den empiriske egenværdi-fordeling konvergere mod den teoretiske
fordeling.
Figur 22: Øverst vises tætheden af egenværdier for korrelationsmatricen
for NYSE77 (lodrette søjler) sammenholdt med Marcenko-Pastur fordelingen (støjbåndet) for Q = 9334/1594 = 5.9 (blå kurve), der ville være
den gældende fordeling hvis aktierne i NYSE77-markedet var ukorrelerede. Det lille indlejrede vindue viser histogrammet i fuld skala. Det ses at
den første egenværdi skiller sig markant ud fra støjbåndet med en værdi
på 340, dvs. ca. 170 gange større end støjbåndets maksimumværdi. Den
2. og 3. egenværdi følger efter med værdier omkring 40. I midten og nederst vises det tilsvarende for hhv. SP2000 og SP84. I alle tre tilfælde ses,
at den empiriske fordeling er væsentlig forskelligt fra Marcenko-Pastur,
hvilket fortolkes på den måde, at markedet indeholder signifikante korrelationer, i modsætning til et komplet ukorreleret marked som det ses på
Figur 21.
Vi kan nu gå skridtet videre og sammenligne Marcenko-Pastur fordelingen
med de empiriske fordelinger for de tre udvalgte markeder, NYSE77, SP84 og
SP2000. Dette er gjort på Figur 22, hvor man ser, at de empiriske egenværdifordelinger adskiller sig radikalt fra Marcenko-Pastur fordelingen. Det kan
dels ses på formen men også (som noteret på figurerne) på cut-off faktoren
k der angiver hvor mange egenværdier, der ligger til højre for λmax og på
korrelationsfaktoren, der til formålet er defineret som forholdet mellem den
højeste empiriske egenværdi λ1 og λmax fra formel 17. Eksempelvis er λ1 = 340
for p
NYSE77 mens λmax = 1.997 som det kan udregnes fra formlen 1 + 1/Q +
2 ∗ (1/Q) hvor Q = 9334/1594.
Vi kan altså konkludere, at NYSE77 har 65 principal-aktiver (svarende til
k = 65) der ligger uden for støj-båndet (Marcenko-Pastur fordelingen) og den
første af disse principal-aktiver har en varians, der er 170 gange højere (angivet på figuren med korrelationsfaktor=170) end man vil finde i et fuldstændig
ukorreleret marked. Med andre ord, der findes i markedet 65 porteføljer, svarende til de første 65 principal-aktiver, som indkapsler al signifikant korrelation
i markedet.
Hermed har vi argumenteret for, at RMT kan anvendes til at give et bud på
en fornuftig cut-off faktor k ud fra betragtningen om, at man kun vil operere
med de principal-aktiver som kan tillægges statistisk signifikans.
56
11 Forecasting
Livet forstås baglæns, men må leves forlæns. - Søren Kierkegaard
Søren Kierkegaards velkendte citat har stor sandhed i relation til risikovurdering. Den sande risiko, den risiko vi virkelig gerne vil kende til, er risikoen for
de hændelser som adskiller sig så meget fra tidligere hændelser, at vi ikke har
en chance for at forudse dem.
Risk managere i små og store virksomheder har ikke den luksus, at kunne
læne sig tilbage sammen med Kierkegaard og se det som et grundvilkår kun
at kunne forstå fortiden. Vurdering af risiko handler i sidste ende altid om
vurdering af fremtidig risiko.
Indtil nu har vi koncentreret os om historiske risikomål, eller udtrykt mere
præcist, om risikomål udregnet for en tidligere periode med kendte markedsdata. Forecasting handler om det modsatte, nemlig om beregning af risikomål
for en fremtidig periode.
For at sætte den rigtige ramme for den følgende diskussion, starter vi med at
tage et kort kig på de metoder og principper, der anvendes i Danske Capital til
at lave forecasting af risiko. Med udgangspunkt heri vil vi give nogle udvalgte
eksempler på hvordan principal-aktiver kan anvendes.
11.1 Forecasting i Danske Capital - et resume
Forecasting i Danske Capital foregår efter en model, der er beskrevet i et
internt notat som er venligt udlånt til denne opgave og endvidere opsummeret
i en artikel i Finans Invest, jf. [9].
Modellen tager som udgangspunkt følgende observationer:
• Afkastfordelinger er ikke normalfordelte. Der er tykke haler og skævhed.
• Ingen egenskaber ved markedet er stationære, alting er i bevægelse. Det
gælder i særdeleshed korrelationer.
• Markedet har regimer, dvs. perioder med forskellige karakteristika. I
det interne notat anvendes en model med fire regimer kaldet expansion,
slowdown, recession og recovery. I artiklen [9] anvendes to regimer kaldet
stresset og ustresset.
• Markedet udviser temporal korrelation, dvs. korrelation mellem afkastet
af et givet aktiv til et tidspunkt og afkastet af samme aktiv til et senere
tidspunkt. Endvidere observeres det, at den temporale korrelation kommer fra de ekstreme hændelser. Hvis man fjerner alle ekstreme hændelser
fra datasættet forsvinder den temporale korrelation næsten.
57
Figur 23: Illustrerer idéen om en regime-model alá den der arbejdes med
i Danske Capital. Man tænker på markedet som et tilstandsrum med et
antal tilstande kaldet regimer, her vises 4 regimer. For hver periode der
går (fx en måned eller et kvartal) er der en given sandsynlighed for at
der sker en overgang til et andet regime.
• De ekstreme hændelser findes med stor overvægt i perioder (dvs. regimer) med recession. Det medfører at regimer med recession vil udvise
temporal korrelation i større grad mens det modsatte vil gøre sig gældende i andre regimer.
• Mange har forsøgt at modellere og parametrisere de multivariate fordelinger i markedet, hvilket har givet anledning til såkaldte copula-teorier.
Da det dels baserer sig på avancerede matematiske koncepter og dermed
bliver svært tilgængeligt for mange investorer og dels giver modeller med
mange frihedsgrader (og dermed stor vilkårlighed i modellen) har Danske Capital valgt en tilgang hvor der i stedet trækkes stikprøver fra de
historiske data. Hermed trækker man altså stikprøver fra den multivariate fordeling som markedet reelt har udvist uden at man behøver en
underliggende matematisk model for fordelingen.
Regime-modellen arbejder med et antal forskellige regimer. I det interne notat
anvendes Recession, Recovery, Expansion og Slowdown. I artiklen [9] anvendes
to regimer kaldet Stresset og Ustresset. Dette kobles med en model for overgangssandsynligheder imellem de forskellige regimer. Denne idé er illustreret
på Figur 23.
For at udregne et givet risikomål for en given portefølje, en given tidshorisont
og en antagelse om et givet intielt regime, bruges regime-modellen først til
at bestemme hvilke regimer der skal betragtes med hvilke vægte over hvilke
58
perioder. Hernæst trækkes historiske hændelser (dvs. afkast) fra de pågældende regimer for den givne portefølje. Hermed får man et bud på en fremtidig
tidsserie for porteføljen, som den vil se ud givet regime-modellen og de historiske data. Slutteligt kan der udregnes alle ønskelige risikomål for porteføljen,
herunder naturligvis varians, VaR og ESF.
Modellen kan udnyttes videre til at lave modeller for optimal allokering af
porteføljer. Princippet er her at man for forskellige risiko-værdier (den ønskede
risiko-profil) laver en afsøgning af markedet som numerisk bestemmer den
portefølje der optimerer afkastet givet risikoen.
11.2 Forecasting med principal-aktiver
Modellen omkring principal-aktiver ændrer ikke ved de grundlæggende diskussioner og problemstillinger man har, når man skal opstille en model for
vurdering af (fremtidig) risiko. Den grundlæggende anke er, at historien ikke
gentager sig selv.
Til gengæld giver principal-aktiverne os en række nye værktøjer og teknikker
til at forstå og analysere markedsdynamikkerne. Dette afsnit vil give nogle
udvalgte eksempler på dette.
Hvis man genlæser afsnittet om forecasting i Danske Capital vil man finde
følgende vigtige analyser:
• Identifikation af regimer.
• Identifikation af hvilket regime man pt. befinder sig i.
• Bestemmelse af overgangssandsynligheder mellem regimer.
• Udregning af afkast for en given portefølje.
• Udregning af realiseret risiko (ex.ante) for en given portefølje for en
historisk periode.
• Estimering af forventet fremtidig risiko (ex.post) for en given portefølje
for et givet regime.
• Bestemmelse af den efficiente rand, dvs. for en givet ønsket risiko (investors ’budget’) find den portefølje der giver højst mulige afkast.
Vi vil her fokusere på spørgsmålet omkring identifikation af regimer. Første
spørgsmål er hvorledes regimer kan bestemmes. Dette spørgsmål giver i sig
selv stof nok til en længere analyse, så vi kommer til at nøjes med en kort
diskussion af emnet. Helt overordnet set må formålet med regimerne være at
finde perioder, hvor markedets karakteristika hvad angår afkast og risiko er
nogenlunde stationære.
59
Lad os her repetere fra afsnit 3 at vi tænker på markedet, som værende beskrevet af n stokastiske variable A1 , . . . , An svarende til de n aktiver i markedet.
Teoretisk set vil A1 , . . . , An være beskrevet af en multivariat fordeling for aktivernes afkast. Vi kan nu formulere kravet ved at sige, at regimer skal dække
perioder hvor den multivariate fordeling er nogenlunde konstant.
Nu er multivariate fordelinger forholdsvis svære at begribe og visualisere. En
særlig gren indenfor statistik kaldet copula-teori har til formål at studere
sådanne fordelinger, blandt andet med henblik på at kunne klassificere og
sammenligne forskellige fordelinger, jævnfør afsnit 5 i [2] og afsnit 5 i [3]. Det
vil føre for vidt at give en formel behandling af spørgsmålet her.
I stedet vil vi simplificere sagen ved at sige, at regimer skal dække perioder
hvor de almindelige korrelationer (dvs. kovarians) er konstante. Korrelationerne er som bekendt givet i form af kovariansmatricen (eller korrelationsmatricen), så man ledes til spørgsmålet, hvordan man herudfra skal identificere
regimerne.
Her kommer principal-aktiverne ind i billedet. Sagen er jo at principal-aktiverne
og de tilhørende egenværdier giver en meget præcis beskrivelse af korrelationerne i markedet. Det første principal-aktiv q1 giver en retning i markedet
for den maksimale varians. En ændring i q1 kan fortolkes som en korrektion i
markedet, hvor korrelationerne har ændret sig, ikke fordi volatiliteten nødvendigvis har ændret sig, men fordi den aktiv-sammensætning (dvs. portefølje)
som giver den største volatilitet (varians) har ændret sig.
Tilsvarende angiver den første egenværdi størrelsen på den maksimale volatilitet i markedet. En ændring i den første egenværdi angiver således en ændring
i markedets volatilitet og dermed et muligt regimeskifte.
Med udgangspunkt i disse ideer kan vi nu lave en model for identificering af
regimer vha. principal-aktiverne som illustreret på Figur 24. Den øverste figur
viser udviklingen af de 6 største egenværdier henover en periode på 30 år. Der
er brugt et glidende vindue på 4 år der skubbes måned for måned. Fx. skal
de første punkter på grafens venstre side hørende til årstallet 1988 fortolkes
således, at der i perioden 1984-1988 har været en markedsdynamik, hvor den
maksimale kovarians (dvs. den maksimale egenværdi) lå et sted mellem 0.01
og 0.02. Det ses også at den højeste egenværdi falder drastisk omkring år 1992
hvorefter den stiger jævnt til ca. år 2004, hvorefter den falder igen.
På den midterste figur er vist et såkaldt egenvektor-delta. Grafen er lavet ved
måned for måned at finde projektionen af den første egenvektor set i forhold
til den første egenvektor måneden før.
projektion til tiden t = q1 (t) · q1 (t − 1)0
Dvs. grafen giver et indtryk af hvilke korrektioner eller chok der sker i markedet. Så længe der er ro i markedet og sålænge den første egenvektor ikke
rigtig ændrer sig henover tid, da ligger projektionen tæt på 1 (projektionen
60
Figur 24: Giver to forskellige overblik over udviklingen i principalaktiverne henover en periode på 30 år. Den øverste graf viser udviklingen
i de 6 første egenværdier. Den midterste graf viser hvilke forskydninger
der har været i det første principal-aktiv fra måned til måned.
er 1 hvis og kun hvis de to vektorer er 100% ens). Når der sker en korrektion
falder projektionen til et sted mellem 0 og 1. På figuren ses således nogle tydelige korrektioner omkring år 1993, igen i år 2003 samt i år 2007 og 2009.
Hvis man sammenligner med den forrige graf vil man se at disse korrektioner
passer ret fint med ændringer i markedsvolatiliteten.
Hermed afrunder vi vores diskussion af principal-aktiver og forecasting. Det
er blevet demonstreret hvorledes modellen åbner op for nye analyse-teknikker
og måder at anskue markedsdynamikkerne på.
61
12 Konklusioner og resultater
Der er givet en definition og gennemført udledning af alle relevante egenskaber ved principal-aktiverne ud fra en central matematisk hovedsætning spektralsætningen. Der er ført beviser for at de er ukorrelerede, at de optimerer variansen og at de udgør en komplet basis for repræsentation af porteføljer
i markedet (se afsnit 4.2 og hovedresultatet (8)). Det er endvidere vist, hvorledes man finder porteføljer på den efficiente rand, når man arbejder med
principal-aktiver (se formel (9)).
Det er vist at man kan anvende principal-aktiver til at lave meget effektive
approximationer af både varians/standardafvigelse (se formel (10)), value-atrisk (se formel (14)) og expected shortfall (se formel (16)).
Der er lavet analyser af historiske data for tre forskellige tidsserier (SP2000,
SP1984 og NYSE77) som strækker sig over et stort antal aktier (op til knap
1600 aktier for NYSE77) og en forholdsvis lang periode - op til 37 år for NYSE77 (se afsnit 3.2). Herunder er der udviklet en række Matlab-programmer
til formålet, hvoraf et uddrag er vist i appendix.
En reduktion i beregningstiden på 99% er generelt muligt for diversificerede
porteføljer i alle tre markeder samtidig med at man nemt opnår en præcision
på omkring 1% for både varians (se afsnit 6.2.3) og expected shortfall (se
afsnit 9.1). Approximation af VaR fungerer lidt dårligere med typiske relative
fejl på omkring 4-5% (se afsnit 8.2) .
Der er gjort rede for en konkret forecasting model, som anvendes i Danske
Capital (se afsnit 11.1). Med udgangspunkt i denne model er der givet et
eksempel på, hvorledes principal-aktiver kan bruges til at udvikle og evt. forbedre denne type forecasting-modeller (se afsnit 11.2).
13 Perspektivering
Vi har med denne opgave kradset i overfladen og taget udvalgte spadestik i
den emnekreds, der åbner sig med introduktionen af principal-aktiver. Den
grundlæggende idé har fra starten været, at finde en ny måde at faktorisere
finansielle markeder på med udgangspunkt i matematisk teori, som allerede
har vist sit værd i mange sammenhænge indenfor finansiering, statistik, fysik
og andre videnskaber.
Det er min overbevisning, at principal-aktiver som forståelsesramme udgør et
vigtigt alternativ til traditionelle faktormodeller. Mange interessante egenskaber er blevet udledt og diskuteret men der er naturligvis mange flere aspekter
som kan og bør undersøges for at få det fulde udbytte af modellen, herunder
kan nævnes følgende:
62
• Det kan undersøges, hvorvidt udledningen af approximationen for VaR
og ESF i afsnit 8 og 9.1 kan underbygges af mere stringente argumenter.
Fx er det forholdsvis oplagt at en multivariat elliptisk fordeling jævnfør
defintion 3.26 i [2] vil have den egenskab, at der er et fast forhold mellem std.afvigelse og fraktil-baserede mål såsom value-at-risk og expected
shortfall. Dermed kan udledningen af (12) forbedres rent teoretisk - uden
at det dog ændrer på approximationens gyldighed og anvendelighed i
virkelighedens verden.
• Det kan undersøges hvor stort et datagrundlag man skal bruge for at
udregne principal-aktiverne op til et givet cut-off k med statistisk signifikans kontra hvor stort et datagrundlag man skal bruge for at udregne
kovarians-matricen med statistisk signifikans. Det skal gerne kunne vises
at de første k principal-aktiver (dvs. både egenværdier og egenvektorer)
kræver et væsentlig mindre datagrundlag end kovariansmatricen som
helhed. Dvs. også her kan der opnås en væsentlig besparelse hvad angår
beregningstid og beregningskompleksitet.
• Som det forklares i [6] og [7] afsnit V.B. er det muligt at rense kovariansmatricen for den støj, som svarer til principal-aktiver i støjbåndet, jævnfør teorien i afsnit 10. Hermed opnår man ifølge disse forfattere bedre
resultater hvad angår allokering af optimale porteføljer.
• Der kan laves reelle målinger der konkret efterviser de besparelser i beregningstid på 99% eller mere, som principelt er påvist i afsnit 6.2.3
• Der kan laves konkrete anvendelser og eftervisninger af formel 9 til fremfinding af den efficiente rand. Det kan eftervises at man kommer frem til
samme efficiente rand som en traditionel Markovitz optimering. Endvidere ville det være interessant at undersøge, om der kan laves effektive
approximationer af den efficiente rand ved at bruge hale-argumenter svarende til dem der blev givet for udledning af varians-approximationen i
afsnit 6.
• Det kan undersøges om man kan føre ideen om principal-aktiver endnu
videre ved at finde de aktiv-sammensætninger i markedet som optimerer ikke variansen men derimod fx VaR−α eller ESFα for et givet α.
På denne måde kunne man måske skabe en endnu bedre basis for at
analysere fraktil-baserede risikomål, jævnfør afsnit VI.C. i [7].
• Man kan gentage mange af de undersøgelser, som er gennemført i denne
opgave, blot hvor kovarians-matricen udregnes som et såkaldt løbende
eksponentielt dæmpet gennemsnit.
Med denne liste af idéer til videre behandling afrundes denne opgave.
63
Litteratur
[1] I. T. Jolliffe,
Principal Component Analysis,
Springer,
Second Edition,
2002.
[2] McNeil, Frey, Embrechts,
Quantitative Risk Management,
Princeton University Press,
First Edition,
2005.
[3] Michael B. Miller,
Mathematics & Statistics for Financial Risk Management,
Wiley,
Second Edition,
2014.
[4] Richard A. Johnson
Miller & Freunds Probability and Statstics for Engineers,
PHI,
Eight Edition.
[5] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
Numerical Recipes in C,
Cambridge,
Second Edition.
[6] Laurent Laloux, Pierre Cizeau, Marc Potters, Jean-Philippe Bouchaud
Random matrix theory and financial correlations ,
Science & Finance,
https://www.math.nyu.edu/faculty/avellane/LalouxPCA.pdf,
World Scientific Publishing Company
[7] Jean-Philippe Bouchaud, Marc Potters
Financial Applications of Random Matrix Theory: a short review,
Science & Finance,
http://arxiv.org/pdf/0910.1205.pdf
[8] Yazann Romahi and Katherine S. Santiago,
Diversification - still the only free lunch?,
https://am.jpmorgan.com/blobcontent/800/973/1383169203651_11_566.pdf,
JPMorgan Asset Management
64
[9] Preben Bertelsen, Rasmus Majborn, Per Søgaard-Andersen,
Risikomåling for langsigtede investorer,
Finans Invest 1/2010,
jf. http://www.finansinvest.dk/Default.aspx?Page=172
[10] Analysegruppen, Aarhus School of Business,
Portfolie Optimization and Matrix Calculation in Excel,
http://nikolaschou.dk/principalaktiver/Portfolio-Optimization.pdf
[11] Wikipedia authors,
Eigenvalue algorithm,
http://en.wikipedia.org/wiki/Eigenvalue_algorithm
[12] Wikipedia authors,
Stock Market Index,
http://en.wikipedia.org/wiki/Stock_market_index
[13] Wikipedia authors,
Risk Parity,
http://en.wikipedia.org/wiki/Risk_parity
[14] Wikipedia authors,
Principal component analysis,
http://en.wikipedia.org/wiki/Principal_component_analysis
Appendices
A Programkode
For at give et indtryk af, hvorledes analyserne i denne opgave er udarbejdet,
vises her to centrale programfiler fra de Matlab-programmer, der er udviklet.
classdef Rates
%UNTITLED4 Encapsulates various utility services around fetching,
%loading and saving rates to desk.
%
properties
end
methods(Static)
function dates = loadTradingDaysChron(period)
%Loads trading days chronologically
dates = flip(Rates.loadTradingDays(period));
end
function dates = loadTradingDays(period)
%%Function loadTradingDays
%Loads trading days for each period between given start and end,
%including both ends. Dates are returned in reverse order (newest
%first)
startYear=period.startYear;
startMonth = period.startMonth;
endYear = period.endYear;
65
endMonth = period.endMonth;
fprintf([’Fetching dates from ’ Utils.formatDate(startYear,startMonth) ’ to ’ Utils.formatDate(endYear, endMonth)]);
dateIntervals = getDateIntervals(startYear, startMonth, endYear, endMonth);
dates = [];
for j=numel(dateIntervals):-1:1
dateObj1 = dateIntervals{j};
ds = Rates.loadTradingDaysForMonth(dateObj1.year, dateObj1.month);
% now concatenate the full time series for a single sec along
% dimension 1
dates = cat(1,dates,ds);
fprintf([’.’]);
end
disp(’ - DONE!’);
end
function cacheFolder = getCacheFolder()
cacheFolder = ’/Users/nikolaschou-pedersen/Documents/MATLAB/cache/’;
end
function n = noTradingDaysMonth(startYear, startMonth)
n = numel(Rates.loadTradingDaysForMonth(startYear, startMonth));
end
function n = noTradingDays(period)
dateIntervals = getDateIntervals(period.startYear, period.startMonth, period.endYear, period.endMonth);
n = 0;
for j=1:numel(dateIntervals)
dateObj1 = dateIntervals{j};
n = n+ Rates.noTradingDaysMonth(dateObj1.year, dateObj1.month);
end
end
function p = getPeriod(y1, m1, y2, m2)
%If only two arguments are provided these will be taken
%startYear and endYear respectively.
if nargin==2
y2=m1;
m1=1;
m2=1;
end
p=struct(’startYear’,y1,’startMonth’,m1,’endYear’,y2,’endMonth’,m2);
end
%%Function loadAll
%% Loads the trading days corresponding to a given period as a column vector.
function dates = loadTradingDaysForMonth(startYear, startMonth)
persistent map;
if isempty(map)
map = containers.Map;
end
key = Utils.formatDate(startYear,startMonth);
if ~map.isKey(key)
folder = [Rates.getCacheFolder() ’_dates/’];
if ~exist(folder,’dir’)
mkdir(folder)
end
dateObj1 = Rates.getDateObj(startYear, startMonth);
dateObj2 = Rates.getDateObj(startYear, startMonth+1);
fileName = [folder,’dates-’,dateObj1.title,’.mat’];
if exist(fileName,’file’)
loadObj = load(fileName);
data = loadObj.data;
else
d1 = datenum(dateObj1.year,dateObj1.month,1);
d2 = datenum(dateObj2.year,dateObj2.month,1);
% Subtract one day so date intervals are non-overlapping
d2 = addtodate(d2,-1,’day’);
%Use IBM to fetch the trading days
data = fetch(Rates.getYahoo(),’IBM’,’Close’, d1, d2);
% only keep date information
data = data(:,1);
save(fileName,’data’);
end
map(key)=data;
end
dates = map(key);
end
function conn = getYahoo()
conn = yahoo(’http://download.finance.yahoo.com’);
end
function rates = loadOrFetchGeneric(conn,sec, period, dataType)
%Loads adjusted prices for the given security
%If the data has been saved to a file, load from there instead
%timer = tic();
folder = Rates.getFolder(sec,dataType);
if ~exist(folder,’dir’)
mkdir(folder)
end
secNameCleaned=strrep(sec,’/’,’_’);
secNameCleaned=strrep(secNameCleaned,’^’,’_’);
fileName = [folder,secNameCleaned,’-’,strrep(Utils.formatPeriod(period),’ to ’,’_’),’.mat’];
66
if exist(fileName,’file’)
temp=load(fileName);
rates=temp.data;
fprintf([’1’ sec]);
else
d1 = Utils.yearMonth2date(period.startYear, period.startMonth);
d2 = Utils.yearMonth2date(period.endYear, period.endMonth);
try
data=fetch(conn,sec,dataType,d1,d2);
fprintf([’2’ sec]);
data=flip(data,1);%Make it chronological
catch exc
data=[];
fprintf([’!’ sec]);
end
save(fileName,’data’);
rates=data;
end
%timer2=toc(timer);
%display(timer2);
end
function rates = fetchAdjusted(sec,period)
data = fetch(yahoo,sec,’Adj Close’,Utils.startDate(period),Utils.endDate(period));
rates = flip(data,1);
end
function title = makeTitle(o)
title = ’’;
if isfield(o,’title’)
title = o.title;
end
if isfield(o,’name’) && isfield(o,’period’)
title=[title o.name ’ ’ Utils.formatPeriod(o.period)];
end
%make recursion
if isfield(o,’parentObject’)
title = [title ’ / ’ Rates.makeTitle(o.parentObject)];
end
end
function variance = getGainsVariance(rateObj)
if ~isfield(rateObj,’gainsCleanedVariance’)
rateObj.gainsCleanedVariance = var(rateObj.gainsChronCleaned);
end
variance = rateObj.gainsCleanedVariance;
end
function plotRates(rateObj,securities,type)
%Utility function to plot rates
if nargin<3
type=’rates’;
end
if isstr(securities)
%Wrapp it in a cell structure
securities = {securities};
end
if iscell(securities)
indexes = Utils.findInCellMany(rateObj.securitiesCleaned,securities);
else
%Assume it is given by a vector of index-positions in the cleaned set
indexes=securities;
securities=rateObj.securitiesCleaned(indexes);
end
if strcmp(type,’rates’)
rates=rateObj.ratesChronCleaned(:,indexes);
elseif strcmp(type,’gains’)
rates=rateObj.gainsChronCleaned(:,indexes);
elseif strcmp(type,’loggains’)
rates=rateObj.loggainsChronCleaned(:,indexes);
else
error(’Unkown type’);
end
period=rateObj.period;
sec=strjoin(securities,’, ’);%[num2str(numel(rateObj.securities)) ’ securities’];
dates=rateObj.dates(1:size(rates,1));
if ~isempty(dates)
plot(dates,rates);
datetick(’x’,’yyyy-mm’);
else
plot(rates);
end
title([’Plot of ’ type ’ (’ sec ’) ’ Utils.formatPeriod(period)]);
zoomAdaptiveDateTicks(’on’);
end
function market = getMarketRates(period)
%Loads rates for the market as a whole based on the S&P index
market= Rates.loadAllAdjusted({’^GSPC’},period);
end
67
function rateObj = loadAllAdjusted(securities, period)
%Loads all rates for all securities for the given period.
%A rateObj is returned
conn = Rates.getYahoo();
timer=tic();
r={};
r.securities=securities;
r.securitiesFailed=containers.Map;
r.securitiesConstrained=containers.Map;%Contains securities not having valid rates for the whole period
display([’Loading ’ num2str(numel(securities)) ’ securities’]);
%Will just use IBM to fetch dates
d1 = Utils.startDate(period);
d2 = Utils.endDate(period);
ibmDates=flip(fetch(conn,’IBM’,’Close’,d1,d2),1);%Use IBM to get the trading days
ibmDates=ibmDates(:,1);
r.noTradingDays=numel(ibmDates);
r.dates=ibmDates;
r.period=period;
result = zeros(r.noTradingDays,numel(securities));
%
for i=1:numel(securities)
fprintf([’ ’ num2str(i) ’/’ num2str(numel(securities)) ’ ’]);
sec = securities{i};
data=Rates.loadOrFetchGeneric(conn,sec,period,’Adj Close’);
if isempty(data)
%Just ignore this security fixing the price at 1,
%think of it like an asset with no movements
result(:,i)=1;
r.securitiesFailed(sec)=i;
else
dataRates=data(:,2);
dataDates=data(:,1);
firstDate=dataDates(1);
firstIndex=find(ibmDates==firstDate,1,’first’);
% *******************************************************
********** MERGING TRADING DAYS ************
% *******************************************************
%As we have no guarantee that ibm and other stocks are
%traded at the same days (for instance, I found that IBM wan
%not traded 27-Sep-1985 whereas AAN was.
%To come around this we will move values day by day
counterIBM=firstIndex;
for j=1:numel(dataDates)
if counterIBM>numel(ibmDates)
%This happens in rare cases where dataDates has
%one more date than ibmDates
break;
end
myDate=dataDates(j);
ibmDate=ibmDates(counterIBM);
if myDate<ibmDate
%We found data for a day not valid for IBM, just
%ignore it and wait for next loop iteration
temp=1;%just to allow debugging
elseif myDate>ibmDate
%IBM has more trading days than the loaded data
%Insert yesterdays rate to keep the rate steady
%Make a loop to increment counterIBM until ibmDate
%catches up with myDate
while(myDate>=ibmDate && counterIBM<numel(ibmDates))
result(counterIBM,i)=result(counterIBM-1,i);%Will never happen in the first iteration so this is safe
counterIBM = counterIBM+1;
ibmDate=ibmDates(counterIBM);
end
else
%Dates must match, copy rate
result(counterIBM,i)=dataRates(j);
counterIBM = counterIBM+1;
end
end
counterIBM = counterIBM-1;%we want it to reflect the last index assigned to
if counterIBM<numel(ibmDates)
%If stock has stopped trading let it have constant rate (to
%avoid high volatility caused by a sudden drop to 0).
result(counterIBM+1:end,i)=result(counterIBM,i);
end
if firstIndex>1
%If stock has not been traded from the beginning let it have
%constant rate until it is starting to trade to avoid high
%volatility caused by a sudden start of trading
result(1:firstIndex-1,i)=dataRates(1);
end
%r.securitiesFailed(sec)=true();
%display([’Fetched end date: ’, datestr(d0), ’, expected end date: ’, datestr(d1)]);
end
end
conn.close();
r.ratesChron=result;
68
r.timer=toc(timer);
rateObj=r;
display(’ ’);
display([’Done - running time ’ num2str(r.timer)]);
end
function ratesObj = loadAll(listOfSecurities, period, options)
%Load all rates for the given period. Rates are returned in
%reverse order, that is, newest first.
%This function will build up a matrix as follows (using 4
%securities as an example):
%-------month 1-----------timerVal = tic;
startYear=period.startYear;
startMonth = period.startMonth;
endYear = period.endYear;
endMonth = period.endMonth;
includedSecurities=containers.Map;
excludedSecurities=containers.Map;
dateIntervals = getDateIntervals(startYear, startMonth, endYear, endMonth);
% We will allocate enough storage and then crop in the end
timeSeriesLength=Rates.noTradingDays(period);
resultTotal = zeros(timeSeriesLength,numel(listOfSecurities));
%Used to keep track of the accumulated number of trading days
%in building up the total matrix
tradingDayCounter = 0;
for j=numel(dateIntervals):-1:1
dateObj1 = dateIntervals{j};
noSecs=size(resultTotal,2);
noTradingDaysThisMonth = Rates.noTradingDaysMonth(dateObj1.year,dateObj1.month);
resultForInterval = zeros(noTradingDaysThisMonth, noSecs);
disp(’ ’);
fprintf([dateObj1.title ’: ’]);
secCounter=1;
for i = 1:numel(listOfSecurities)
secu = listOfSecurities{i};
if ~excludedSecurities.isKey(secu)
data = Rates.loadOrFetchRatesForMonth(secu, dateObj1.year, dateObj1.month, options.buildArray);
if (options.buildArray)
if Rates.isValidRates(data) && numel(data)==noTradingDaysThisMonth
resultForInterval(:,secCounter) = data;
secCounter = secCounter + 1;
includedSecurities(secu)=true();
else
% Remove this security from the matrix
excludedSecurities(secu)=true();
resultForInterval(:,secCounter) = [];
resultTotal(:,secCounter) = [];
end
end
end
end
% now concatenate the full time series for a single sec along
% dimension 1
if (options.buildArray)
resultTotal(tradingDayCounter+1:tradingDayCounter+noTradingDaysThisMonth,:) = resultForInterval;
tradingDayCounter = tradingDayCounter + noTradingDaysThisMonth;
end
end
rates = resultTotal;
disp(’ ’);
ratesObj = {};
ratesObj.ratesChron = flip(rates,1);
ratesObj.rates = rates;
ratesObj.excludedSecurities = excludedSecurities;
remove(includedSecurities,keys(excludedSecurities));
ratesObj.includedSecurities = includedSecurities;
ratesObj.period = period;
ratesObj.elapsedTime = toc(timerVal);
ratesObj.dateIntervals = dateIntervals;
end
function [ rates ] = loadOrFetchRatesForMonth( secu, startYear, startMonth, returnData )
%loadOrFetchSecurity Loads or fetches the rates for a given security for
%one month
%
First it is checked whether the security can be found in the file
%
system. Otherwise it is fetched and then saved to desk.
%If returnData is false and a data file is found in the cache, this
%function will return quickly.
folder = Rates.getFolder(secu);
if ~exist(folder,’dir’)
mkdir(folder)
end
dateObj1 = Rates.getDateObj(startYear, startMonth);
dateObj2 = Rates.getDateObj(startYear, startMonth+1);
secName = strrep(secu,’/’,’_’);
fileName = [folder,secName,’-’,dateObj1.title,’.mat’];
if exist(fileName,’file’)
if (returnData)
loadObj = load(fileName);
69
data = loadObj.data;
else
% Ok, no need to do more
data = [];
end
fprintf([’ 1’ secu]);
else
d1 = datenum(dateObj1.year,dateObj1.month,1);
d2 = datenum(dateObj2.year,dateObj2.month,1);
% Subtract one day so date intervals are non-overlapping
d2 = addtodate(d2,-1,’day’);
try
data = fetch(Rates.getYahoo(),secu,’Close’, d1, d2);
% throw away date information
data = data(:,2);
% fprintf([’Fetched from yahoo: ’,secu,’ ’,dateObj1.title]);
fprintf([’ 2’ secu]);
catch exc
data = [];
fprintf([’ !’ secu]);
% fprintf([’Data cannot be loaded for security ’, secu, ’. Used zeros instead.’]);
end
save(fileName,’data’);
end
rates = data;
end
function val=getDateObj(startYear,startMonth)
d = datenum(startYear,startMonth,1);
val = {};
vec =datevec(d);
val.date = d;
val.year = vec(1);
val.month = vec(2);
paddedMonth = sprintf(’%02d’, val.month);
mytitle = [num2str(val.year),’-’,paddedMonth];
val.title=mytitle;
end
function val=isValidRates(data)
%at least 5 trading days
if numel(data)>5 && data(1)>0 && data(end)>0
val=true();
else
val=false();
end
end
function folder = getFolder(sec,type)
%Returns the folder containing cached data for the given
%security sec
if nargin==1
type=’Close’;
end
folder = [Rates.getCacheFolder(), ’_’,strrep(type,’ ’,’’),’/’,sec, ’/’];
end
function fileName = getMetafile(sec)
fileName = [Rates.getFolder(sec) ’meta.mat’];
end
function meta = getMeta(sec)
%Loads or inits a meta object for the given security
metafile = Rates.getMetafile(sec);
if ~exist(metafile,’file’)
meta = {};
else
loadObj = load(metafile);
meta = loadObj.meta;
end
end
function saveMeta(sec,meta)
%Saves the meta object for the given security
metafile = Rates.getMetafile(sec);
save(metafile,’meta’);
end
function yearMonth = findFirstValidMonth(sec)
%Determines which month is the first valid month for which rates
%across the whole month can be fetched for the given security
%The ’first valid month’ M is characterized by the following facts:
%1. Non-zero rates are found for M
%2. A number of non-zero rates are found corresponding to the
%number of banking days for that month
%3. Bullet 1 and 2 are not satisfied for the previous month.
meta = Rates.getMeta(sec);
if ~isfield(meta,’firstValidMonth’)
meta.firstValidMonth = Rates.findFirstValidMonthUtil(sec);
Rates.saveMeta(sec,meta);
70
end
yearMonth = meta.firstValidMonth;
end
function yearMonth = findFirstValidMonthUtil(sec)
%Finds the first valid month using a kind of log2-algorithm
persistent allMonths;
if isempty(allMonths)
%apparently yahoo only delivers data since 1963
allMonths = getDateIntervals(1963,1,2014,1);
end
N = numel(allMonths);
index = round(N/2);
leftEnd = 1;
rightEnd = N;
counter = 0;
if Rates.isValidMonth(sec,allMonths{1}.year, allMonths{1}.month)
%First check if it has been defined from the very beginning
index = 1;
else
while(true())
%display([’Investigate index ’ num2str(index) ]);
m1 = allMonths{index};
isValid = Rates.isValidMonth(sec,m1.year, m1.month);
if isValid
m0 = allMonths{index - 1};
if ~Rates.isValidMonth(sec, m0.year, m0.month)
%if the previous month is invalid we have found the
%right month
break;
end
end
%If we reach this point we know that we did not find the
%first valid month
if isValid
%go left
rightEnd = index;
increment = round((index - leftEnd)/2);
index = index - increment;
if index <= 1
index = 1;
break;
end
else
%go right
leftEnd = index;
increment = round((rightEnd - index)/2);
index = index + increment;
if index >= N
index = N;
break;
end
end
if counter > N
%Just to avoid infinite loops
error(’Could not find the first valid month’);
end
counter = counter + 1;
end
end
yearMonth = allMonths{index};
end
function b = isValidMonth(sec, y, m)
rates = Rates.loadOrFetchRatesForMonth(sec, y, m, true);
if min(rates) == 0
%A month with a rate of 0 cannot be valid
b = false();
else
noTradingDays = Rates.noTradingDaysMonth(y,m);
if numel(rates) == noTradingDays
b = true();
else
b = false();
end
end
end
function [ result ] = filterSecurities( secus, period)
%Filters a list of securites by the criteria that rates should
%be available for the given yearMonth.
result = cell(1,numel(secus));
counter = 1;
for i=1:numel(secus)
sec = secus{i};
%First check whether the data is actually valid on the
%given yearMonth
isStartValid = Rates.isValidMonth(sec, period.startYear, period.startMonth);
isEndValid = Rates.isValidMonth(sec, period.endYear, period.endMonth);
if isStartValid && isEndValid
result{counter} = secus{i};
71
counter = counter + 1;
end
end
if counter<numel(secus)
% shrink it
result(counter:end)=[];
end
end
end
end
classdef Analysis
%UNTITLED2 Summary of this class goes here
%
Detailed explanation goes here
properties
end
methods(Static)
function rateObj=prepare(rateObj)
%Make some further preparation of a rateObj needed for the analysis
rateObj.dates=Rates.loadTradingDaysChron(rateObj.period);
rateObj.datesForPlot=Utils.formatDatesForPlot(rateObj.dates);
rateObj.securities=keys(rateObj.includedSecurities);
rateObj.gainsChron = Analysis.findGains(rateObj.ratesChron);
rateObj.diffsChron = rateObj.ratesChron(2:end,:) - rateObj.ratesChron(1:end-1,:);
end
function rateObj=prepareAdjusted(rateObj, troubledSecurityNames)
%troubledSecurityIndexes: An array of security indexes
%troubledSecurityNames
%Make some further preparation of a rateObj needed for the analysis
rateObj.datesForPlot=Utils.formatDatesForPlot(rateObj.dates);
%Now clean up the rates matrix
indexes = find(prod(rateObj.ratesChron)==0);%Indexes of columns having at least one zero (which should never happen)
for i=1:numel(troubledSecurityNames)
secIndex=Utils.findInCell(rateObj.securities,troubledSecurityNames{i});
if ~isempty(secIndex)
indexes(numel(indexes)+1) = secIndex;
end
end
indexes = unique(indexes);
rateObj.ratesChronCleaned=rateObj.ratesChron;
rateObj.ratesChronCleaned(:,indexes)=[];
rateObj.securitiesCleaned=rateObj.securities;
rateObj.securitiesCleaned(indexes)=[];
rateObj.securityIndexesTroubled=indexes;
rateObj.gainsChronCleaned = Analysis.findGains(rateObj.ratesChronCleaned);
rateObj.loggainsChronCleaned = log(rateObj.gainsChronCleaned+1);
rateObj.diffsChronCleaned = rateObj.ratesChronCleaned(2:end,:) - rateObj.ratesChronCleaned(1:end-1,:);
rateObj.n=numel(rateObj.securitiesCleaned);
rateObj.T=size(rateObj.dates,1);
end
function rateObj=prepareUni(uni,period,includedSecuritiesList)
%Takes the result of buildUniverse yahoo command and aligns it
%with the structure of other rateObjs
rateObj={};
rateObj.period=period;
rateObj.securities=includedSecuritiesList;
rateObj.ratesChron=uni(:,2:end);
rateObj.dates=uni(:,1);
rateObj.datesForPlot=Utils.formatDatesForPlot(rateObj.dates);
rateObj.gainsChron = Analysis.findGains(rateObj.ratesChron);
rateObj.diffsChron = rateObj.ratesChron(2:end,:) - rateObj.ratesChron(1:end-1,:);
end
function showEigen(corrMatrix,period,figNr,figTitle, securities)
%Plots the eigenvalues of a given corrMatrix or covMatrix for a
%given period
figure(figNr);
[eigvec,eigval]=Utils.eig2(corrMatrix);
noBins = round(numel(eigval)/3);
eigval=eigval(1:10);
hist(eigval,noBins);
clear title xlabel ylabel;
title([’Eigenvalues of ’ figTitle ’ ’ Utils.formatPeriod(period)]);
figure(figNr+1);
plot(sort(eigvec(:,1)));
title([’First eigenvector of ’, figTitle, ’ ’, Utils.formatPeriod(period)]);
[max0,index0]=max(abs(eigvec(:,1)));
display([’First eigenvector of ’, figTitle, ’ dominated by ’, securities{index0}]);
end
function full(rates, period, securities)
%%Makes an analysis of different
gains = Analysis.findGains(rates);
72
%loggains = log(gains);
%Analysis.showEigen(corr(rates),period,1,’rates’);
Analysis.showEigen(corr(gains),period,1,’correlation of gains’, securities);
Analysis.showEigen(cov(gains),period,3,’covariance of gains’, securities);
%Analysis.showEigen(corr(loggains),period,5,’log of gains’);
%Analysis.showEigen(corr(exp(gains)),period,7,’exp of gains’);
end
function drawRandomMatrix(sizeT,sizeN)
r2=rand(sizeT,sizeN);
eigvalrand=eig(corr(r2));
figure(9);
hist(eigvalrand,20);
title([’Eigenvalues of correlation of random matrix of size ’ num2str(sizeT) ’X’ num2str(sizeN)]);
end
function gains = findGains(rates)
%Determines the gains from a timeseries of rates
gains = rates(2:end,:)./rates(1:end-1,:) - 1;
end
function analysis = riskForecasting(rateObj, quantile, portfolios, pcaPeriod, figNo)
%
analysis={};
analysis.pcaAnalysis = Analysis.PCA(rateObj,pcaPeriod,figNo);
analysis.periodLength=pcaPeriod.endYear-pcaPeriod.startYear;
forecastingLength = rateObj.period.endYear - pcaPeriod.endYear;
analysis.portfolioValueAtRiskErrorMean = zeros(forecastingLength+1,1);
analysis.portfolioValueAtRiskErrorStdDev = zeros(forecastingLength+1,1);
analysis.portfolioESFErrorMean = zeros(forecastingLength+1,1);
analysis.portfolioESFErrorStdDev = zeros(forecastingLength+1,1);
%Slide through the years
for i = 0:forecastingLength
analysisTemp=Analysis.valueAtRisk(analysis.pcaAnalysis,quantile,portfolios,Rates.getPeriod(pcaPeriod.startYear+i,pcaPeriod.endYear+i));
analysis.portfolioValueAtRiskErrorMean(i+1)=analysisTemp.portfolioValueAtRiskErrorMean;
analysis.portfolioValueAtRiskErrorStdDev(i+1)=analysisTemp.portfolioValueAtRiskErrorStdDev;
analysis.portfolioESFErrorMean(i+1)=analysisTemp.portfolioESFErrorMean;
analysis.portfolioESFErrorStdDev(i+1)=analysisTemp.portfolioESFErrorStdDev;
end
end
function analysis = valueAtRisk(pcaAnalysis, quantile, portfolios, forecast, cutoff)
%
Finds the historical value at risk measure for a given portfolio
%fractile: A number close to 0, e.g. 0.05 for the 5% quantile
%analysis: Expects an analysis object as returned from the PCA
%
function
%portfolios: A n x m matrix where each row corresponds to a vector of stock weights summing to 1.
%nrYearsToForecast: Can either be a number measuring the number
%
of years to forecast ahead of the period used for PCA or it
%
can be a period.
%init
if nargin<3
portfolios=Utils.createPortfolios(numel(pcaAnalysis.eigenValues),20);
end
if nargin<4
forecast = 2;
end
if ~isfield(forecast,’startYear’)
%Assume it is an integer indicating the number of years to forecast
p = pcaAnalysis.periodUsedForPCA;
forecast = Rates.getPeriod(p.endYear,p.endMonth,p.endYear + forecast,p.endMonth);
end
if nargin<5
cutoff = 10;
end
analysis={};
analysis.portfolios=portfolios;
analysis.quantile=quantile;
analysis.cutoff = cutoff;
analysis.periodForecasted=forecast;
[analysis.minIndex, analysis.maxIndex] = Utils.findDateInterval(pcaAnalysis.rateObj.dates,forecast);
min=analysis.minIndex;
max=analysis.maxIndex;
analysis.pcaAnalysis=pcaAnalysis;
analysis.parentObject=pcaAnalysis;
analysis.rateObj=pcaAnalysis.rateObj;
%find VaR of assets
gains = pcaAnalysis.rateObj.gainsChronCleaned(min:max,:);
[analysis.valueAtRisk, analysis.ESF] = Analysis.valueAtRiskUtil(gains,quantile);
analysis.gainsVariance = var(gains);
stdDev = sqrt(analysis.gainsVariance);
analysis.valueAtRiskFromVariance = stdDev * norminv(quantile,0,1);
%VaR and ESF for eigenvectors should be found from the
%PCA-based period
[analysis.eigenValueAtRisk, analysis.eigenValueESF] = Analysis.valueAtRiskUtil(pcaAnalysis.eigenGains,quantile);
73
%analysis.eigenValueAtRisk has one row with VaR calculated
%historically for each eigenvector
analysis.eigenValueAtRiskFromVariance = sqrt(pcaAnalysis.eigenValues) * norminv(quantile,0,1);
if Utils.equalPeriods(forecast,analysis.pcaAnalysis.periodUsedForPCA)
analysis.title=[’ \alpha =’ num2str(analysis.quantile*100) ’%’];
else
analysis.title=[’Forecast for ’ Utils.formatPeriod(forecast) ’, \alpha =’ num2str(analysis.quantile*100) ’%’];
end
analysis.portfolioGains = gains*transpose(portfolios);%Gains for the forecast period projected on each portfolio
[analysis.portfolioValueAtRisk,analysis.portfolioESF] = Analysis.valueAtRiskUtil(analysis.portfolioGains,quantile);
analysis.eigenProjections = portfolios * pcaAnalysis.eigenVectors;%Each row will contain projection constants for the respective portfolio
%At this point the top row of eigenProjections is:
%|alpha1 alpha2 ... alpha_n|
%being projections of
%portfolios(1,:), i.e. the ’first’ portfolio, on the first
%eigenVector.
%Exploit that eigenvectors are uncorrelated, hence
%VaR(portfolio) = sqrt(sum of (eigenvalueVaR^2 x projection_weight^2))
analysis.portfolioValueAtRiskComposedFromEigenVectors = sqrt(analysis.eigenValueAtRisk.^2 * (analysis.eigenProjections.^2)’);
%Now calculate with cutoff
analysis.valueAtRiskTail = mean(analysis.eigenValueAtRisk(cutoff+1:end).^2);
tailAddition = analysis.valueAtRiskTail*sum(analysis.eigenProjections(:,cutoff+1:end).^2,2); %column vector with entries for each portfolio
analysis.portfolioValueAtRiskComposedFromEigenVectorsWithCut = sqrt( tailAddition’ + analysis.eigenValueAtRisk(1:cutoff).^2 * (analysis.eigenPro
analysis.portfolioESFComposedFromEigenVectors = sqrt(analysis.eigenValueESF.^2 * (analysis.eigenProjections.^2)’);
%Now calculate with cutoff
analysis.ESFTail = mean(analysis.eigenValueESF(cutoff+1:end).^2);
tailAddition2 = analysis.ESFTail*sum(analysis.eigenProjections(:,cutoff+1:end).^2,2); %column vector with entries for each portfolio
analysis.portfolioESFComposedFromEigenVectorsWithCut = sqrt( tailAddition2’ + analysis.eigenValueESF(1:cutoff).^2 * (analysis.eigenProjections(:
analysis.portfolioValueAtRiskError = (analysis.portfolioValueAtRiskComposedFromEigenVectors-analysis.portfolioValueAtRisk) ./ analysis.portfolio
analysis.portfolioValueAtRiskErrorMean = mean(analysis.portfolioValueAtRiskError);
analysis.portfolioValueAtRiskErrorStdDev = sqrt(var(analysis.portfolioValueAtRiskError));
analysis.portfolioValueAtRiskErrorWithCut = (analysis.portfolioValueAtRiskComposedFromEigenVectorsWithCut-analysis.portfolioValueAtRisk) ./ anal
analysis.portfolioValueAtRiskErrorMeanWithCut = mean(analysis.portfolioValueAtRiskErrorWithCut);
analysis.portfolioValueAtRiskErrorStdDevWithCut = sqrt(var(analysis.portfolioValueAtRiskErrorWithCut));
analysis.portfolioESFError = (analysis.portfolioESFComposedFromEigenVectors-analysis.portfolioESF) ./ analysis.portfolioESF;
analysis.portfolioESFErrorMean = mean(analysis.portfolioESFError);
analysis.portfolioESFErrorStdDev = sqrt(var(analysis.portfolioESFError));
analysis.portfolioESFErrorWithCut = (analysis.portfolioESFComposedFromEigenVectorsWithCut-analysis.portfolioESF) ./ analysis.portfolioESF;
analysis.portfolioESFErrorMeanWithCut = mean(analysis.portfolioESFErrorWithCut);
analysis.portfolioESFErrorStdDevWithCut = sqrt(var(analysis.portfolioESFErrorWithCut));
%Analysis.plotValueAtRisk(analysis, figNo);
end
function plotValueAtRisk(analysis, figNo)
%%Plot value at risk for original assets
figure(figNo);
a=analysis.valueAtRiskFromVariance;
b=analysis.valueAtRisk;
temp = cat(1,a,b);
temp = (transpose(temp));
hlines=plot(temp);
mytitle=[’Value-at-risk (’ Rates.makeTitle(analysis) ’)’];
title(mytitle);
xlabel(’Aktienr’);
ylabel(’VaR’);
set(hlines(1),’Displayname’,’Udregnet fra standardafvigelser’);
set(hlines(2),’Displayname’,’Udregnet via historisk distribution’);
legend(’Location’,’north’);
%%Plot value at risk for eigen vectors
figure(figNo+1);
a=analysis.eigenValueAtRiskFromVariance;
b=analysis.eigenValueESF;
temp = cat(1,a,b);
temp = (transpose(temp));
hlines=plot(temp);
mytitle=[’Expected Shortfall for Eigenvectors (’ Rates.makeTitle(analysis) ’)’];
title(mytitle);
xlabel(’Eigenvector’);
ylabel(’ESF’);
set(hlines(1),’Displayname’,’Udregnet fra standardafvigelser’);
set(hlines(2),’Displayname’,’Udregnet via historisk distribution’);
legend(’Location’,’north’);
%%Plot value at risk for eigen vectors
figure(figNo+2);
a=analysis.eigenValueAtRiskFromVariance;
74
b=analysis.eigenValueAtRisk;
temp = cat(1,a,b);
temp = (transpose(temp));
hlines=plot(temp);
mytitle=[’Value-at-risk for Eigenvectors (’ Rates.makeTitle(analysis) ’)’];
title(mytitle);
xlabel(’Eigenvector’);
ylabel(’VaR’);
set(hlines(1),’Displayname’,’Udregnet fra standardafvigelser’);
set(hlines(2),’Displayname’,’Udregnet via historisk distribution’);
legend(’Location’,’north’);
if size(analysis.portfolios,1)>1
%%Plot value at risk for portfolios
figure(figNo+3);
a=analysis.portfolioValueAtRisk;
b= analysis.portfolioValueAtRiskComposedFromEigenVectors;
temp = Utils.easyCat(a,b);
hlines=plot(temp);
mytitle=[’Portfolio value-at-risk (’ Rates.makeTitle(analysis) ’)’];
title(mytitle);
xlabel(’Portfolio number’);
ylabel(’VaR’);
set(hlines(1),’Displayname’,’Udregnet fra historiske data’);
set(hlines(2),’Displayname’,’Udregnet via PCA sum’);
legend(’Location’,’north’);
%%Plot ESF for portfolios
figure(figNo+4);
a=analysis.portfolioESF;
b= analysis.portfolioESFComposedFromEigenVectors;
temp = Utils.easyCat(a,b);
hlines=plot(temp);
mytitle=[’Portfolio Expected Shortfall (’ Rates.makeTitle(analysis) ’)’];
title(mytitle);
xlabel(’Portfolio number’);
ylabel(’ESF’);
set(hlines(1),’Displayname’,’Udregnet fra historiske data’);
set(hlines(2),’Displayname’,’Udregnet via PCA sum’);
legend(’Location’,’north’);
%%Plot VaR delta between historical VaR and PCA based VaR
if size(analysis.portfolios,1)>50
figure(figNo+5);
hist(analysis.portfolioValueAtRiskError, numel(analysis.portfolioValueAtRiskError)/30);
mytitle=[’Portfolie VaR relativ fejl, gns=’ Utils.formatPct(analysis.portfolioValueAtRiskErrorMean) ’, std.afv.=’ Utils.formatPct(analys
title(mytitle);
xlabel(’Relativ fejl’);
ylabel(’Hyppighed’);
Utils.formatAxisPct(’x’);
figure(figNo+6);
hist(analysis.portfolioESFError, numel(analysis.portfolioESFError)/30);
mytitle=[’Portfolie ESF relativ fejl, gns=’ Utils.formatPct(analysis.portfolioESFErrorMean) ’, std.afv.=’ Utils.formatPct(analysis.portf
title(mytitle);
xlabel(’Relativ fejl’);
ylabel(’Hyppighed’);
Utils.formatAxisPct(’x’);
figure(figNo+7);
hist(analysis.portfolioValueAtRiskErrorWithCut, numel(analysis.portfolioValueAtRiskErrorWithCut)/30);
mytitle=[’Portfolie VaR relativ fejl, cutoff=’ num2str(analysis.cutoff) ’, gns=’ Utils.formatPct(analysis.portfolioValueAtRiskErrorMeanW
title(mytitle);
xlabel(’Relativ fejl’);
ylabel(’Hyppighed’);
Utils.formatAxisPct(’x’);
figure(figNo+8);
hist(analysis.portfolioESFErrorWithCut, numel(analysis.portfolioESFErrorWithCut)/30);
mytitle=[’Portfolie ESF relativ fejl, cutoff=’ num2str(analysis.cutoff) ’, gns=’ Utils.formatPct(analysis.portfolioESFErrorMeanWithCut)
title(mytitle);
xlabel(’Relativ fejl’);
ylabel(’Hyppighed’);
Utils.formatAxisPct(’x’);
end
end
end
function [valueAtRisk, ESF] = valueAtRiskUtil(gains, quantile)
%quantile: Typically 0.95. Measures the VaR quantile
n=size(gains,1);
m=size(gains,2);
quantileIndex=round(n*(1 - quantile));
valueAtRisk = zeros(1,m);
ESF = zeros(1,m);
for i=1:m
75
s=sort(gains(:,i));
valueAtRisk(1,i)=-s(quantileIndex);
ESF(1,i) = -mean(s(1:quantileIndex));
end
end
function analysis = PCA(rateObj, periodUsedForPCA, figNo)
%This analysis will inspect the market evolution of the portfolios
%corresponding to the first few principal components of the
%whole period considered.
if nargin<=1
periodUsedForPCA=rateObj.period;
end
if nargin<=2
figNo = 1;
end
d1 = Utils.startDate(periodUsedForPCA);
d2 = Utils.endDate(periodUsedForPCA);
[minIndex, maxIndex] = Utils.findDateInterval(rateObj.dates, periodUsedForPCA);
analysis={};
analysis.rateObj=rateObj;
analysis.parentObject=rateObj;
analysis.periodUsedForPCA=periodUsedForPCA;
analysis.minIndex=minIndex;
analysis.maxIndex=maxIndex;
analysis.gainsUsedForPca=rateObj.gainsChronCleaned(minIndex:maxIndex,:);
[V,E]=Utils.eig2(cov(analysis.gainsUsedForPca));
analysis.eigenValues=E;
analysis.eigenVectors=V;
temp=sum(V);
n = numel(temp);
%Normalize so that the eigenvectors have sum = 1
%which means we can think of it as portfolio weights
analysis.eigenVectorsNormalized = V ./ (ones(n,n)*diag(temp));
analysis.eigenGains = analysis.gainsUsedForPca * analysis.eigenVectors;%These are the gains for the PCA-period projected on the eigen vectors
analysis.eigenGainsFullPeriod = rateObj.gainsChronCleaned * analysis.eigenVectors;%These are the gains for the whole period projected on the eige
analysis.eigenRatesFullPeriod=rateObj.ratesChronCleaned * analysis.eigenVectors;
%Makes no sense : analysis.eigenRatesNormalized=rateObj.ratesChronCleaned * analysis.eigenVectorsNormalized;
%Makes no sense: analysis.eigenGainsNormalized = analysis.gainsUsedForPca * analysis.eigenVectorsNormalized;
%Makes no sense : analysis.eigenGainsCalculated = Analysis.findGains(analysis.eigenRates);
market = Rates.getMarketRates(rateObj.period);
firstEigen=analysis.eigenRatesFullPeriod(:,1);
ratio=mean(market.ratesChron ./ firstEigen);
analysis.marketRates = market.ratesChron/ratio;
analysis.title=[’PCA on ’ Utils.formatPeriod(analysis.periodUsedForPCA)];
end
function pcaPlot(analysis, figNo)
try
close(figNo);
catch
end
rateObj=analysis.rateObj;
figure(figNo);
ratesCombined=cat(2, analysis.marketRates, analysis.eigenRatesFullPeriod(:,1:3));
hlines = plot(analysis.rateObj.dates,ratesCombined);
datetick(’x’);
mytitle=[’Prisudvikling for 1. og 2. egenvektor - ’ Utils.formatPeriod(rateObj.period) ’ baseret på data fra ’ Utils.formatPeriod(analysis.period
title(mytitle);
xlabel(’Tid’);
ylabel(’Pris’);
set(hlines(1),’Displayname’,’S&P index’);
set(hlines(2),’Displayname’,’1. egenvektor’);
set(hlines(3),’Displayname’,’2. egenvektor’);
set(hlines(4),’Displayname’,’3. egenvektor’);
legend(’Location’,’north’);
%%
figure(figNo+1);
hlines=semilogy(Utils.easyCat(var(analysis.eigenGainsFullPeriod),analysis.eigenValues));
title([’Varians af afkast af egenaktiverne for hele perioden vs. egenværdier ’ Rates.makeTitle(analysis.rateObj)]);
xlabel(’Egenaktiv nummer’);
ylabel(’Afkast’);
set(hlines(1),’Displayname’,’Historisk varians af egenaktiver’);
set(hlines(2),’Displayname’,’Egenværdier’);
legend(’Location’,’north’);
end
function analysis=varianceApproximation(pcaAnalysis, portfolioWeights, cutoff)
%pcaAnalysis: The result of calling PCA()
%portfolioWeights: A portfolio represented by a row-vector of weights, that
%is the vector elements must sum to 1.
%It is assumed that the portfolio is continously rebalanced in
%order for the weigths to be valid throughout the period
%cutoff: A number or vector of numbers controlling the number of
%eigenvectors used to approximate the variance. If a number
%between 0 and 1 (excluding 1) is used, it is considered as a
76
%percentage of the total number. If the number is 1 or above it
%is consideres as an absolute number. If it is a vector of
%numbers the analysis will be made for each cutoff.
%if size(portfolioWeights,1)>1
% error(’This analysis is designed for a single portfolio’);
%end
analysis={};
analysis.portfolioGains=pcaAnalysis.rateObj.gainsChronCleaned*transpose(portfolioWeights);
analysis.portfolios = portfolioWeights;
analysis.varianceHistorical=var(analysis.portfolioGains);
analysis.pcaAnalysis=pcaAnalysis;
analysis.parentObject=analysis.pcaAnalysis;
analysis.rateObj=pcaAnalysis.rateObj;
lambda = diag(pcaAnalysis.eigenValues);
analysis.projections=portfolioWeights*pcaAnalysis.eigenVectors; %beta factors, 1xn vector
analysis.varianceByPca=analysis.projections * lambda * analysis.projections’; %variance by pca using all egenvectors
N = numel(pcaAnalysis.eigenValues);
analysis.cutoff=Utils.either(numel(cutoff)==1 && cutoff(1)<1, round(numel(pcaAnalysis.eigenValues)*cutoff), cutoff);
nrCutoffs = numel(analysis.cutoff);
analysis.varianceByPcaApprox = zeros(size(portfolioWeights,1),nrCutoffs);
sumOfStockVariance = sum(Rates.getGainsVariance(pcaAnalysis.rateObj));
for i=1:nrCutoffs
K = analysis.cutoff(i);
sumOfInitialEigenValuesProjections = (analysis.projections(:,1:K).^2 * pcaAnalysis.eigenValues(1:K)’);
if K<N
fraction = 1/(numel(pcaAnalysis.eigenValues)-K);
sumOfInitialEigenValues = sum(pcaAnalysis.eigenValues(1:K));
lambdaTail = fraction * (sumOfStockVariance - sumOfInitialEigenValues);
vectorLengthOfPortfolioWeights = sum(portfolioWeights.*portfolioWeights,2);
sumOfInitialBetaSquare = sum(analysis.projections(:,1:K) .* analysis.projections(:,1:K),2);
betaTailSquare = vectorLengthOfPortfolioWeights - sumOfInitialBetaSquare;
else
lambdaTail=0;
betaTailSquare=0;
end
%
alphaCut = analysis.projections .* ((1:N)<=K); %alpha cut off
analysis.lambdaTail = mean(pcaAnalysis.eigenValues(K+1:end));
w = portfolioWeights;
analysis.varianceByPcaApprox(:,i) = sumOfInitialEigenValuesProjections + lambdaTail*betaTailSquare;
%
%
end
histVar = repmat(analysis.varianceHistorical’,1,nrCutoffs);
analysis.varianceByPcaApproxError=abs(analysis.varianceByPcaApprox./histVar - 1);
analysis.stdDevApproxError=abs(sqrt(analysis.varianceByPcaApprox./histVar) - 1);
%Average approximation error across all portfolios and all k’s
analysis.stdDevApproxErrorMean = mean(mean(analysis.stdDevApproxError));
%Average (std.dev. of approximation error across portfolios for fixed k) across k
analysis.stdDevApproxErrorMeanStdDev = mean(std(analysis.stdDevApproxError));
analysis.stdDevApproxCutoff=[];
analysis.stdDevApproxCutoffTreshold=[];
%Find the least cutoff K0 such that no approximation error is greater
%than some limit for any k>=K0
for i=1:8
treshold=10*10^(-i);
analysis.stdDevApproxCutoff(i) = Analysis.findTreshold(analysis.stdDevApproxError,treshold);
analysis.stdDevApproxCutoffTreshold(i) = treshold;
end
end
function result = findTreshold(data, limit)
result = find(max(data)>limit,1,’last’)+1;
if isempty(result)
result = 0;
end
end
function varianceAnalysis=varianceApproximationTotal(rateObj,cutoff, portfolios, portfolioName)
pca=Analysis.PCA(rateObj,rateObj.period,1);
%Analysis.pcaPlot(pca1984,4);
%variance1984=Analysis.varianceApproximation(pca1984,portfolios,4);
if nargin<3
portfolios = Utils.createPortfolios(numel(rateObj.securitiesCleaned),1);
end
varianceAnalysis=Analysis.varianceApproximation(pca,portfolios,cutoff);
Analysis.varianceApproximationPlot(varianceAnalysis,1,portfolioName);
Utils.saveFigures(1:3, [’Variance-approx-’ rateObj.name ’-’ portfolioName]);
end
function totalAnalysis = varianceApproximationTotal2(rateObj, indexes)
a={};
if nargin<2
a.indexes = cat(2,1:30,32:2:60, 65:5:rateObj.n, rateObj.n);
77
else
a.indexes=indexes;
end
portfolio1 = diag(ones(1,rateObj.n));
portfolio4 = Utils.createCircularPortfolios(rateObj.n,4);
portfolio8 = Utils.createCircularPortfolios(rateObj.n,8);
portfolio16 = Utils.createCircularPortfolios(rateObj.n,16);
portfolio32 = Utils.createCircularPortfolios(rateObj.n,32);
portfolio64 = Utils.createCircularPortfolios(rateObj.n,64);
%portfolio(1,1)=1;
portfolioDiversified = Utils.createPortfolios(rateObj.n,100);
a.approx1 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio1 ,’Enkelt-aktier’);
a.approx4 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio4 ,’4 aktier’);
a.approx8 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio8 ,’8 aktier’);
a.approx16 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio16 ,’16 aktier’);
a.approx32 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio32 ,’32 aktier’);
a.approx64 = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolio32 ,’64 aktier’);
a.approxDiversified = Analysis.varianceApproximationTotal(rateObj, a.indexes, portfolioDiversified ,’Diversificeret’);
totalAnalysis=a;
end
function totalAnalysis = varianceApproximationTotal3(varAnals)
fields={’approx1’,’approx4’,’approx8’,’approx16’,’approx32’,’approx64’,’approxDiversified’};
a={};
a.fixedPrecision001=zeros(1,1);
a.fixedPrecision0001=zeros(1,1);
a.fixedPrecision001Pct=zeros(1,1);
a.fixedPrecision0001Pct=zeros(1,1);
for i=1:numel(varAnals)
for j=1:numel(fields)
analysis = varAnals(i);
fieldName = fields(j);
varAnalysis = analysis.(fieldName{1});
nrEigenvalues=varAnalysis.rateObj.n;
a.fixedPrecision001(i,j) = varAnalysis.stdDevApproxCutoff(3);
a.fixedPrecision0001(i,j) = varAnalysis.stdDevApproxCutoff(4);
a.fixedPrecision001Pct(i,j)= a.fixedPrecision001(i,j)/nrEigenvalues;
a.fixedPrecision0001Pct(i,j)= a.fixedPrecision0001(i,j)/nrEigenvalues;
end
end
totalAnalysis = a;
end
function varianceApproximationPlot(analysis, figNo, portfolioName)
nrPortfolios=size(analysis.portfolios,1);
mytitle=[’PCA approximation af sigma / ’ portfolioName Utils.either(nrPortfolios>1,[’ / #porteføljer=’ num2str(nrPortfolios)] ,’’) ’ (’ Rates.ma
if nrPortfolios == 1
figure(figNo);
result = Utils.easyCat(ones(1,size(analysis.varianceByPcaApprox,2))*sqrt(analysis.varianceHistorical),sqrt(analysis.varianceByPcaApprox));
hlines=plot(analysis.cutoff,result);
title(mytitle);
xlabel(’Antal egenværdier medregnet’);
ylabel(’Standardafvigelse’);
set(hlines(1),’Displayname’,’Historisk std.afv.’);
set(hlines(2),’Displayname’,’Std.afv. udregnet ved egenværdi-approximation’);
Utils.showLegend();
end
figure(figNo+1);
hlines=plot(analysis.cutoff,analysis.stdDevApproxError);
title(mytitle);
xlabel(’Antal egenværdier medregnet’);
ylabel(’Standardafvigelse fejl’);
set(hlines(1),’Displayname’,’Fejl ved egenværdi-approximation’);
Utils.formatAxisPct(’y’);
if nrPortfolios > 1
figure(figNo+2);
s=analysis.stdDevApproxError;
result=cat(1,mean(s),mean(s)+2*sqrt(var(s)),min(s),max(s));
hlines = plot(analysis.cutoff,result’);
title(mytitle);
xlabel(’Antal egenværdier medregnet’);
ylabel(’Standardafvigelse fejl’);
set(hlines(1),’Displayname’,’Fejl ved egenværdi-approximation’);
Utils.formatAxisPct(’y’);
set(hlines(1),’Displayname’,’Gns. fejl’);
set(hlines(2),’Displayname’,’Gns. fejl 95% konfidens’);
set(hlines(3),’Displayname’,’Minimum fejl’);
set(hlines(4),’Displayname’,’Maximum fejl’);
Utils.showLegend();
end
end
function analysis = eigenEvolution(rateObj, periodLenInYears, figNo, useCovariance, useType)
%%
78
% --- Find the evolution in the top eigenvalues using a sliding 6-month span
% Loop over all months, 13 years each having 12 months
analysis = {};
if strcmp(useType,’gains’)
fullmatrix = rateObj.gainsChronCleaned;
elseif strcmp(useType,’diffs’)
fullmatrix = rateObj.diffsChronCleaned;
elseif strcmp(useType,’loggains’)
fullmatrix = rateObj.loggainsChronCleaned;
elseif strcmp(useType,’rates’)
fullmatrix = rateObj.ratesChronCleaned;
else
error([’Uknown type ’ useType]);
end
period = rateObj.period;
%This contains the number of retained eigenvalues
noEigValues = 6;
%The length of the period for which the correlation matrix is built
periodLen = 251*periodLenInYears;
noSecurities = size(fullmatrix,2);
%Determine the counterMax from the periodLength
%counterMax will be the number of samples we make
%Generally we will make a sample each month
counterMax = numel(getDateIntervals(period.startYear,period.startMonth,period.endYear,period.endMonth));
%The sliding window should be subtracted
counterMax = counterMax - periodLenInYears*12;
%Time scale is vertical, eigenvalues are flipped horizontal
eigvalEvolution = zeros(counterMax, noEigValues+1);
%Keep a copy of the first eigenvector throughout the whole
%evolution
eigvecEvolutionFull = zeros(counterMax, noSecurities);
dates=Rates.loadTradingDaysChron(period);
%Timescale will be vertical
dateTicks=zeros(counterMax,1);
%Will contain the first eigenvector for each time span
disp(’**** Now starting the evolution analysis **** ’)
fprintf([’Running up to ’ num2str(counterMax-1) ’: ’]);
for i=0:counterMax-1
fprintf([num2str(i) ’ ’]);
%%
%first day in the period for this loop
startNo = 1+i*21;
%last day in the period for this loop
endNo = startNo + periodLen;
%We use the last date of the period as tick
dateTicks(i+1) = dates(endNo);
slidingData=fullmatrix(startNo:endNo,:);
if useCovariance
cotemp = cov(slidingData);
else
cotemp = corr(slidingData);
end
[eigvec2,eigval2] = Utils.eig2(cotemp);
eigvalEvolution(i+1,2:end) = transpose(eigval2(1:noEigValues));
%eigvalEvolution(i+1,1) = sum(eigval2);
eigvecEvolutionFull(i+1,:) = Analysis.normalizeFirstEigvec(transpose(eigvec2(:,1)));
end
%Now make the analysis of the eigenvector evolution
eigvecEvolutionMean = mean(eigvecEvolutionFull,1);
%Keep track of the projection of the first eigenvector against
%itself in the preceeding period
eigvecEvolutionDelta = zeros(counterMax,1);
%Keep track of the distance from the mean
eigvecEvolutionDiffusion = zeros(counterMax,1);
for i=1:counterMax
eigVec = transpose(eigvecEvolutionFull(i,:));
if i==1
eigvecEvolutionDelta(i,1)=1;
else
eigVecPreceeding = eigvecEvolutionFull(i-1,:);
eigvecEvolutionDelta(i,1)=Analysis.project(eigVecPreceeding,eigVec);
end
eigvecEvolutionDiffusion(i,1) = Analysis.project(eigvecEvolutionMean,eigVec);
end
analysis.eigvecEvolutionDelta=eigvecEvolutionDelta;
analysis.eigvecEvolutionDiffusion=eigvecEvolutionDiffusion;
analysis.eigvalEvolution=eigvalEvolution;
analysis.dates=dates;
analysis.dateTicks = dateTicks;
analysis.rateObj = rateObj;
analysis.eigvecEvolutionFull=eigvecEvolutionFull;
analysis.periodLenInYears = periodLenInYears;
disp(’Done’);
%analysis.title = [Utils.either(useCovariance,’covariance’,’correlation’) ’ of ’ useType ’ ’ Utils.formatPeriod(period) ’ (sliding window=’ num
analysis.title = [Utils.formatPeriod(period) ’ (sliding window=’ num2str(periodLenInYears) ’y)’];
Analysis.eigenEvolutionPlot(analysis, figNo);
end
79
function analysis=eigenEvolutionTotal(rateObj)
analysis = Analysis.eigenEvolution(rateObj,4,1,true,’gains’);
Utils.saveFigures(1:3,[’Evolution-’ analysis.rateObj.name ’-’ num2str(analysis.periodLenInYears) ’y’]);
end
function eigenEvolutionPlot(analysis,figNo)
%%
figure(figNo);
plot(analysis.dateTicks,analysis.eigvalEvolution);
mytitle = [’Eigenvalue evolution of ’ analysis.title];
title(mytitle);
datetick(’x’,’yyyy-mm’);
zoomAdaptiveDateTicks(’on’);
%%
figure(figNo+1);
plot(analysis.dateTicks,transpose(analysis.eigvecEvolutionDelta));
mytitle = [’Eigenvector delta of ’ analysis.title];
title(mytitle);
datetick(’x’,’yyyy-mm’);
zoomAdaptiveDateTicks(’on’);
%%
figure(figNo+2);
plot(analysis.dateTicks,transpose(analysis.eigvecEvolutionDiffusion));
mytitle = [’Eigenvector diffusion of ’ analysis.title];
title(mytitle);
datetick(’x’,’yyyy-mm’);
zoomAdaptiveDateTicks(’on’);
end
function p = project(eigvec1,eigvec2)
%%
%Finds the projection of eigvec1 against eigvec2. If either eigvec1 or eigvec2
%is found to point in the "opposite" direction they are
%reversed to the same direction. This is done by taking the
%absolute value of eigvec1*eigvec2
p = abs(eigvec1*eigvec2);
end
function result = normalizeFirstEigvec(eigvec)
%Normalizes the first eigenvector in the sense that the
%direction of the vector is changed to have dominating positive
%values. If the eigenvector is mixed positive and negative this
%does not make fully sense, but often the first eigenvector in
%a financial market is dominated by either positive or negative
%coordinates, not mixed
noPositive = sum(arrayfun(@(x)x>0,eigvec));
if noPositive<numel(eigvec)/2
result = -1 * eigvec;
else
result = eigvec;
end
end
end
end
80