Predavanje 25: Primene, Glava 7

FlipView[{Text[Style[" Branko Urošević, imqfbrankourosevic@gmail.com ",20,"Text",Italic]],Text[Style[" Faculty of Economics, University of Belgrade, NBS and CESIfo Munich",20,"Text", Italic]]}]

7. Statistička obrada podataka u WL

WL je izuzetno pogodan kako za grafičko predstavljanje podataka (to smo već videli) tako i za njihovo statističko modeliranje i obradu. U ovom poglavlju ćemo ukratko izložiti kako se u WL radi statističko testiranje hipoteza,  linearna regresija, i analiza vremenskih serija.

Ovde je pre svega cilj da se intuitivno nauči kako da primenimo WL na rešavanje takvih problema, bez previše obaziranja na neke od važnih tehničkih detalja kada nam oni nisu krajnje neophodni. Te tehničke detalje možete naučiti na kursevima ekonometrije i vremenskih serija ili naći u knjigama ili na internetu.

1. Statističko testiranje hipoteza

Uvodna razmatranja

Kada imamo neki skup podataka, interesuje nas kakva funkcija raspodele ili model stohastički proces mu najbolje odgovara, kakva je “prava” očekivana vrednost raspodele iz koje je izvučen taj skup podataka, kakva je varijansa, da li postoji statistička veza između nekih grupa podataka itd.

Odgovori na ova pitanja ne mogu se dati sa sigurnošću. Naime, možemo da kažemo sa kakvom verovatnoćom možemo da odbacimo neku hipotezu ali ne i da tvrdimo da je nešto svakako tačno (ovaj poslednji iskaz bio bi deterministički). Različiti testovi mogu da daju različite rezultate. U onoj meri u kojoj veliki broj različitih testova ukazuje na isti zaključak mi smo sigurniji da smo u pravu. Ovde ćemo da vidimo kako se u WL vrše neki statistički testovi. Pitanja na koja ćemo želeti da odgovorimo, između ostalog, su:

Kako da verifikujemo da je srednja vrednost neke grupe podataka jednaka (manja ili veća) od neke zadate vrednosti u statističkom smislu?
Kako da uporedimo srednje vrednosti dve ili više grupe podataka?
Kako testirati kolika je varijansa nekog uzorka?
Kako da proverimo zavisnost između dve ili više grupa podataka?
Kako da proverimo da li je neki uzorak izvučen iz neke zadate raspodele (recimo normalne)?

Testiranje osnovnih hipoteza vezanih za srednju vrednost

Da bismo proverili da li neka grupa podataka ima srednju vrednost koja se statistički ne razlikuje od 0, mi možemo da postavimo nultu hipotezu da je srednja vrednost jednaka 0. Onda je alternativna hipoteza da srednja vrednost nije jednaka 0.

Ono što možemo da izračunamo (na više različitih načina) jeste verovatnoća sa kojom možemo da prihvatimo (ili odbacimo) nultu hipotezu. Stoga, da bismo testirali hipotezu, postoje funkcije koje računaju upravo verovatnoću da treba da prihvatimo osnovnu hipotezu (odnosno da odbacimo alternativnu). Što je veća ta verovatnoća, to smo sigurniji u našem zaključku.

Hipotezu o srednjoj vrednosti uzorka možemo da proverimo putem LocationTest funkcije.

Razmotrimo uzorak od 1000 slučajno generisanih brojeva iz standardne normalne raspodele.

uzorak=RandomVariate[NormalDistribution[],1000];

Uzoračka sredina i histogram izgledaju ovako:

Mean[uzorak]

SmoothHistogram[uzorak]

Hajde da proverimo da li se srednja vrednost statistički ne razlikuje od nule. To radimo pomoću

LocationTest[uzorak,0]

Testiranje uzorka iz standardne normalne raspodele

To znači da je po nekom standardnom testu (po defoltu) ovaj uzorak napravljen iz raspodele koja ima srednju vrednost 0 sa verovatnoćom od (šta god da smo dobili tu) posto. U pitanju je p-vrednost o čemu ćemo nešto više reći kasnije. Što je p-vrednost veća to je manja verovatnoća da možemo da odbacimo originalnu hipotezu.

Kada bismo imali verovatnoću ispod 10 posto, morali bismo da odbacimo nultu hipotezu sa pouzdanošću od 90 posto. Možemo da istovremeno obavimo više testova kako bismo sigurniji u našim zaključcima:

LocationTest[uzorak,Automatic,{"TestDataTable",All}]

Ne ulazeći u detalje toga kakvi su ovi testovi vidimo da ne možemo statistički da odbacimo nultu hipotezu da je srednja vrednost uzrorka jednaka nuli.

Testiranje srednje vrednosti uzoraka iz drugačijih normalnih raspodela

Napravimo sada dva uzorka, jedan iz raspodele sa teorijskom srednjom vrednošću 5, a drugi sa teorijskom srednjom vrednošću 1.

uzorak1=RandomVariate[NormalDistribution[5,2],1000];

uzorak2=RandomVariate[NormalDistribution[1,1],1000];

Testirajmo sada nultu hipotezu da je srednja vrednost oba uzorka jednaka 5. To bismo uradili ovako:

LocationTest[#,5]&/@{uzorak1,uzorak2}

Vidimo da za prvi uzorak nultu hipotezu ne možemo da odbacimo dok za drugu možemo da odbacimo sa velikom pouzdanošću (što naravno nije iznenađujuće ako se setimo iz kojih raspodela smo generisali uzorke). Naravno, ukoliko bismo testirali hipotezu da su srednje vrednosti jednake 1, dobili bismo da moramo da odbacimo tu hipotezu za prvi ali ne i za drugi uzorak.

LocationTest[#,1]&/@{uzorak1,uzorak2}

Testiranje za koliko je srednja vrednost jednog uzorka veća od srednje vrednosti drugog

Proverimo sada hipotezu da je srednja vrednost prvog uzorka za 4 veća od srednje vrednosti drugog. To se piše ovako:

LocationTest[{uzorak1,uzorak2},4]

Vidimo da po defoltnom testu ne možemo da odbacimo tu hipotezu. Hajde da proverimo to isto za sve testove odjednom. Uradimo to ovako:

h0=LocationTest[{uzorak1,uzorak2},4,"HypothesisTestData"]

Evo kako možemo da izvučemo iz toga neke važne statistike:

h0["PValue","T"]

h0["TestStatistic","T"]

Možemo da više testova prikažemo i grafički, i evo kako.

resultati=h0["TestData",All];

testovi=h0["AllTests"];

BarChart[resultati[[All,2]],ChartLabels-> Placed[testovi,Center],BarOrigin-> Left]

Testiranje srednje vrednosti kada je hipoteza data preko nejednakosti

Do sada smo testirali samo hipoteze da je srednja vrednost ili razlika srednje vrednosti jednaka nečemu. Umesto toga, možemo da testiramo hipoteze koristeći nejednakosti.

Na primer, kako da testiramo hipotezu da je razlika u srednjim vrednostima između dva uzoraka veća ili jednaka 4?  Tako što ćemo definisati alternativnu hipotezu preko "Less". Naime:

LocationTest[{uzorak1,uzorak2},4,AlternativeHypothesis-> "Less"]

Dakle, ne možemo odbaciti hipotezu da je razlika bar 4.

Interesatno je da možemo da menjamo i nivo statističke pouzdanosti koristeći sledeću opciju

Clear[hipoteza0]

hipoteza0=LocationTest[{uzorak1,uzorak2},4,AlternativeHypothesis-> "Less",SignificanceLevel->0.05]

Rezultat testiranja prikazan rečima

Napravimo ponovo uzorak:

uzorak1=RandomVariate[NormalDistribution[5,2],1000];

LocationTest[uzorak1,5]

Proverimo sada hipotezu da je srednja vrednost jednaka 5 na nivou signifikanosti 0.1, 0.05, 0.01, i 0.00001. Ono što nas zanima jeste da prikažemo rezultat testiranja rečima, a ne brojkom.

{ℋ1,ℋ2,ℋ3,ℋ4}=LocationTest[uzorak1,5,"HypothesisTestData",SignificanceLevel->#]&/@{0.1, 0.05, 0.01, 0.00001};

ℋ1["TestConclusion"]

ℋ2["TestConclusion"]

ℋ3["TestConclusion"]

ℋ4["TestConclusion"]

Nastavak

Nije teško videti da ukoliko testiramo da je srednja vrednost jednaka 3 umesto 5 imali bismo:

{ℋ1,ℋ2,ℋ3,ℋ4}=LocationTest[uzorak1,3,"HypothesisTestData",SignificanceLevel->#]&/@{0.1, 0.05, 0.01, 0.00001};

ℋ1["TestConclusion"]

ℋ2["TestConclusion"]

ℋ3["TestConclusion"]

ℋ4["TestConclusion"]

Testiranje hipoteza vezanih za varijansu

Varijansa je druga veoma važna karakteristika slučajnih uzoraka. Njen nivo se testira vrlo slično kao kod srednje vrednosti. Funkcija za to je VarianceTest.

Pretpostavimo da želimo da generišemo uzorak od 1000 standardnih normalnih promenljivih.

uzorak5=RandomVariate[NormalDistribution[],1000];

Izračunajmo sada uzoračku varijansu:

Variance[uzorak5]

Da li se ona statistički razlikuje od teorijske vrednosti 1? To možemo da testiramo koristeći

VarianceTest[uzorak5,1]

Da li možemo da odbacimo da uzorak koji smo generisali dolazi iz raspodele sa varijansom jednakom 1? To zavisi od uzorka.

Nastavak

Probajmo sada nešto drugo. Da li možemo da odbacimo da je varijansa zapravo 1.1?

VarianceTest[uzorak5,1.05]

VarianceTest[uzorak5,1.1]

A šta je sa 0.9?

VarianceTest[uzorak5,0.9]

Na taj način, možemo naći interval u kome bi razumno morala da bude prava varijansa.

Kako testirati jednakost varijanse dva različita uzorka?

Ništa prostije. Naime, pretpostavimo da smo napravili dva uzorka. Oba imaju istu teorijsku varijansku jednaku 1.

uzorak10=RandomVariate[NormalDistribution[2,1],1000];

uzorak11=RandomVariate[NormalDistribution[3,1],1000];

SmoothHistogram[{uzorak10,uzorak11},PlotStyle->{Thin,Thick}]

Ovako testiramo jednakost varijansi:

VarianceEquivalenceTest[{uzorak10,uzorak11}]

Vidimo da po defolt testu ne možemo da odbacimo da su varijanse identične. Šta više, to ne možemo da uradimo ni za celu seriju testova:

VarianceEquivalenceTest[{uzorak10,uzorak11},{"TestDataTable",All}]

Kako proveriti da li su dve ili više grupa podataka nezavisne jedna od druge

Prilikom rada sa podacima veoma je važno znati da li su podaci statistički nezavisni ili zavisni. Od toga zavisi kakvi modeli imaju smisla a kakvi ne.

U WL postoji funkcija IndependenceTest koja to proverava. Ona radi veoma slično kao i druge funkcije za testiranje hipoteza koje smo već videli. Nulta hipoteza jeste da su grupe podataka nezavisne a alternativna je da nisu.

Razmotrimo Binormalnu raspodelu sa korelacijom jednakom 0. U tom slučaju očekujemo da dobijemo statistički nezavisne podatke za dve kolone. Hajde da to generišemo.

RandomVariate[BinormalDistribution[0],10]

Ukoliko hoćemo da te tačke predstavimo preko dve kolone brojeva, moramo da ovu listu transponujemo. Definišimo dve varijable (kolone) kao:

v1

{v1,v2}=RandomVariate[BinormalDistribution[0],200]\[Transpose];

Možemo da izračunamo korelaciju:

Nastavak

Correlation[{v1,v2}\[Transpose]]

Vidimo da je koeficijent korelacije mali. Da li je to, međutim, statistički različito od nule, odnosno da li su dve kolone v1, v2 nezavisne? Primetimo razliku u sintaksi ovde. Naime, nema vitičastih zagrada kao u prethodnim primerima!

IndependenceTest[v1,v2]

Po defoltu dobijamo da se nezavisnost ne može odbaciti. Hajde da sada probamo da dobijemo odgovore na sve testove:

IndependenceTest[v1,v2,{"TestDataTable",All}]

Očigledno da ni po kom kriterijumu ne možemo da odbacimo nezavisnost ove dve varijable. Primetimo da male vrednosti koeficijenta korelacije govore o tome da, možda, ne možemo da ih razlikujemo statistički od nule.

Naći raspodelu koja najbolje opisuje neku grupu podataka

Razmotimo vrednost indeksa SP500 u 2015. To je indeks koji opisuje kretanje cena 500 najvećih kompanija u SAD i služi kao proxy za tržišni portfolio. U tom slučaju:

sp500=FinancialData["SP500",{{2015}}][[All,2]]

sp500=Last[FinancialData["SP500",{2015}]\[Transpose]]

Ukoliko internet ne radi importujemo fajl sp500indexu2015.dat.

Izračunajmo sada logaritamske prinose

prinossp500=Log[Most[sp500]/Rest[sp500]];

SmoothHistogram[prinossp500,PlotRange->All]

Ovako dobijamo raspodelu koja po određenom kriterijumu predstavlja najbolji fit.

FindDistribution[prinossp500,4]

Nastavak

Možemo sada da vidimo kako izgleda ta procenjena raspodela verovatnoće

PDF[FindDistribution[prinossp500],x]

Plot[PDF[LogisticDistribution[-6.230220651160469`*^-6,0.0054621290928360215`],x],{x,-0.06,0.06},PlotRange-> All]

Testirajmo sada u kojoj meri procenjena raspodela odgovara podacima

h0=DistributionFitTest[prinossp500,NormalDistribution[0.000035069669720310966`,0.010339083477691573`],"HypothesisTestData"]

h0["TestDataTable",All]

NormalDistribution[0.000035069669720310966`,0.010339083477691573`]

h0=DistributionFitTest[prinossp500,FindDistribution[prinossp500],"HypothesisTestData"]

h0["TestDataTable",All]

Vidimo da ta raspodela ne može da odbaci nultu hipotezu, naime da je originalna raspodela iz koje su podaci dobijeni gornja LaplaceDistribution.

2, Linearna regresija u WL

Da li su žene diskriminisane sa platama u SAD?

Pretpostavimo da želimo da ustanovimo da li su žene diskriminisane u odnosu na muškarce što se tiče plate. Jedan prost, ali ne i preterano pouzdan, način da to ustanovimo jeste da napravimo regresiju plate po satu rada u odnosu na varijablu koja je jednaka 0 ukoliko je u pitanju žena i 1 ukoliko je u pitanju muškarac.

Onda probamo da objasnimo platu pomoću modela u kome se kao objašnjavajuća varijabla pojavljaje ta koja određuje da li je neko muškarac ili žena. Ukoliko je ta varijabla signifikantna, onda to može da znači da to što je neko muškarac ili ne utiče na platu. Razume se, iz toga je prerano zaključiti da li postoji diskriminacija. Ali, biće neka indikacija da bar postoji šansa da je tako.

Podatke u WL unosimo, na primer, preko Import funkcije, ovako (koristimo Insert i onda File path:

wagesdata=Import["C:\\Users\\mathematica\\Desktop\\Podaci za WL Primene\\wages1.dat","Table"]\[Transpose];

Pogledajmo dimenzije ove tabele:

Dimensions[wagesdata]

Definisanje varijabli

Svaki red u ovoj velikoj matrici odgovara jednoj promenljivoj. Na primer, treći red je:

wagesdata[[4]]

Definišimo sada varijable za platu po satu rada (wages). To je zavisna promenljiva koju želimo da objasnimo.

wagesdata[[4]]

Odbacimo prvi element (ime varijable) i definišimo:

wages=Rest[wagesdata[[4]]];

Broj godina iskustva (exper), da li je neko muškarac (1) ili žena (0) (to je varijabla male) i najzad broj godina školovanja (school) datu su na sledeći način:

male=Rest[wagesdata[[2]]]

exper=Rest[wagesdata[[1]]];

school=Rest[wagesdata[[3]]];

Linearna regresija sa jednom nezavisnom promenljivom

Probajmo sada da objasnimo platu po satu rada samo na osnovu toga da li je ispitanik muškarac ili žena. Formiramo matricu u kojoj idu prvo objašnjavajuće variable i na kraju variabla koja treba da se objasni:

data1={male,wages}\[Transpose];

Vidimo da prvi podatak odgovara ženi čija je plata 6.31 dolar na sat, drugi podatak odgovara ženi čija je plata 5.47 dolara na sat itd.

Take[data1,30]

Evo sada i kako se fituju ovi podaci koristeći LinearModelFit

model1=LinearModelFit[data1,xmale,xmale]

Output je objekat koji sadrži velike količine korisnih informacija o modelu ukoliko znamo kako da ih koristimo. Prvo, vidimo da, saglasno modelu, koeficijent uz male varijablu ne izgleda jednak nuli.

Ako je significantan, to bi značilo da u srednjem muškarci zarađuju 1.12 dolara više od žena (srednja plata žena je 5.15 otprilike po satu rada, to je konstantni koeficijent). Jednačinu modela možemo da izvučemo ovako:

Normal[model1]

Nastavak o linearnoj regresiji sa jednom nezavisnom promenljivom

Korisno je, takođe, da vidimo kako izgledaju intervali poverenja za procenjene parametre. To radimo ovako:

model1["ParameterConfidenceIntervalTable"]

Primetimo da su oba procenjena koeficijent statistički značajna i pozitivna (intervali poverenja na nivou 95 posto ne sadrže nulu). Druge korisne stvari koje možemo videti iz modela jeste:

model1["ANOVATable"]

Vidimo da koeficijenti nisu ne samo pojedinačno nego ni zajedno jednaki nuli (veoma mala P vrednost odnosno veoma veliki F koeficijent)

Još jedna standardna karakteristika modela jeste (R^2) (udeo varijacije koji objašnjava varijabla male). Ona je niska:

model1["RSquared"]

Linearna regresija sa dve nezavisne promenljive

Naravno, nije teško zapitati se: da li bi rezultat bio signifikantan ukoliko bismo u obzir uzeli, recimo, i godine iskustva? Jer, može da se desi da su muškarci iskusniji na poslu i, stoga, bolje plaćeni? U tom slučaju, ako kao objašnjavajuću variablu dodamo exper može doći do promene.

data2={exper,male,wages}\[Transpose];

Take[data2,30]

Evo kako bi se fitovao ovakav model pomoću LinearModelFit

model2=LinearModelFit[data2,{xexper,xmale},{xexper,xmale}]

model2a=LinearModelFit[data2,{xexper,(xexper^2),xmale},{xexper,xmale}]

Dakle, novi model koji smo dobili jeste:

Normal[model2a]

model2a["ParameterConfidenceIntervalTable"]

model2a["RSquared"]

Normal[model2]

Opet moramo da proverimo statističku značajnost. Za pojedinačne varijable dobijamo sledeće intervale poverenja:

model2["ParameterConfidenceIntervalTable"]

Vidimo da je exper zapravo, pojedinačno gledano, nesignificantna. Sa druge strane, male je i dalje statistički značajna.

model2["ANOVATable"]

Vidimo da kolektivno nisu sve beznačajne. Sa druge strane, dodavanje ove jedne varijable doprinosi malo poboljšanju (R^2).

model2["RSquared"]

Nastavak regresije sa dve objašnjavajuće promenljive

Kada imamo više varijabli, onda je bolje koristiti Adjusted (R^2) (on uzima u obzir da smo možda veštački povećali (R^2) samo zato što smo povećali broj promenljivih)

model2["AdjustedRSquared"]

Probajmo, sada, umesto exper da u model ubacimo godine obrazovanja. Naime,

data3={school,male,wages}\[Transpose];

Evo sada i kako se fituje LinearModelFit

model3=LinearModelFit[data3,{xschool,xmale},{xschool,xmale}]

model3a=LinearModelFit[data3,{xschool,(xschool^2),xmale},{xschool,xmale}]

Normal[model3a]

Normal[model3]

Da li možda više godina školovanja muškaraca objašnjava razliku u plati muškaraca i žena? Ne liči da je tako. Vidimo da su sada svi koeficijenti statistički značajni:

model3a["ParameterConfidenceIntervalTable"]

model3a["AdjustedRSquared"]

model3["ParameterConfidenceIntervalTable"]

Primetimo, takođe, da je Adjusted (R^2) značajno unapređen u odnosu na prethodne modele.

Mean[school]//N

model3["AdjustedRSquared"]

Regresija sa 3 objašnjavajuće promenljive

Da li bi se stvari dodatno popravile ukoliko bismo stavili sve tri objašnjavajuće promenljive?

data4={school,exper,male,wages}\[Transpose];

model4=LinearModelFit[data4,{xschool,xexper,xmale},{xschool,xexper,xmale}];

Normal[model4]

Proverimo sada statističku značajnost varijabli:

model4["ParameterConfidenceIntervalTable"]

Primetimo da je sada i xexper postala značajna! Takođe,

model4["AdjustedRSquared"]

Vidimo da smo dobili model koji ima relativno pristojan Adjusted (R^2).

Analiza reziduala

Da bismo bili sigurni da naša regresiona analiza ima smisla, moramo da proverimo svojstva reziduala. Kako izgledaju reziduali? Jedan način da se proveri pretpostavka o normalnosti grešaka bila bi da se napravi takozvani QuantilePlot.

On crta empirijske kvantile serije grešaka i kvantile normalne raspodele. Ukoliko bi kvantili grešaka ležali na pravoj liniji, oni bi bili normalno raspodeljeni. Očigledno, to ovde nije slučaj.

model4["FitResiduals"]

QuantilePlot[model4["FitResiduals"]]

Nije teško ustanoviti da su greške nisu normalno raspodeljene

h0=DistributionFitTest[model4["FitResiduals"],NormalDistribution[a,b],"HypothesisTestData"]

h0["TestDataTable",All]

Još o rezidualima

Proverimo sada da li su reziduali nezavisni u odnosu na varijable sistema. Vidimo ovako da nema očigledne zavisnosti.

ListPlot[model4["FitResiduals"]]

Ali, proverimo sada nezavisnost reziduala u odnosu na 3 objašnjavajuće varijable:

IndependenceTest[model4["FitResiduals"],#]&/@{school,exper,male}

Vidimo da se na nivou signifikantnosti od 5 posto ili manje ne može odbaciti nezavisnost. Ukoliko, pak, koristimo 10 posto kao kriterijum, onda se nezavisnost može praktično odbaciti. Znači, ovaj model je na granici prihvatljivosti ovakav kakav je.

Model prinosa na hartije od vrednosti: Fama-French faktori

Ideja je da se analizira koji su glavni faktori koji utiču na prinose akcija na tržištu. Na njih, pre svega, utiče kretanje na tržištu u celini, kao i neki drugi faktori rizika koji se zovu Fama-French faktori (dobili su ime po dvojici ekonomista koji su ih uveli u razmatranje). Podatke unosimo iz fajla capm.dat.

Tu imamamo podatke o prinosima na neke portfolije hartije od vrednosti kao i neke faktore rizika koji treba da budu objašnjavajuće varijable. Podaci su mesečni, od januara 1960 do oktobra 2010 (610 meseci). Podaci su, tipično, dati kao prinosi u odnosu na nerizičnu stopu (u procentima): dodatni prinos = prinos- rf. Takođe imamo i podatke o nerizičnoj stopi rf. Još jedna varijabla koju možemo da definišemo jeste januarska promenljiva. Ona je jednaka 1 kada je mesec januar, i 0 u ostalim slučajevima.

capmdata = Transpose[Import["C:\\Users\\mathematica\\Desktop\\Podaci za WL Primene\\capm.dat","Table"]];

capmdata[All,1]

Grid[Map[Style[#,20,Bold]&,{{"Portfolio","Objašnjenje"},{"foodrf","dodatni prinos na indeks industrije hrane"},{"durblrf","dodatni prinos na indeks nepotrošne robe"},{"cnstrrf","dodatni prinos na indeks građevinarstva"},{"rmrf","dodatni prinos na tržišni portfolio"},{"rf","prinos na nerizičnu aktivu"},{"smb","dodatni prinos na male kompanije"},{"hml","dodatni prinos na hml faktor"},{"jan","dummy varijabla za januar"}},{2}],Frame->All]

Take[constrrf,20]

{constrrf,durblrf,foodrf,hml,jan,rf,rmrf ,smb}=Rest[capmdata[[#]]]&/@Range[2,9];

CAPM model

Probajmo da objasnimo dodatne prinose na foodrf indeks (to jest na indeks industrije hrane u Americi). Počnimo od najprostijeg modela, naime CAPM. Dodatni prinos se onda pretpostavlja da zavisi samo od  dodatnog prinosa na tržišni portfolio rmrf. Odredimo  koeficijente regresije (alfa i beta) za ovaj portfolia.

Formirajmo prvo matricu podataka. Pošto su svi podaci u procentima, moramo da ih sve podelimo sve brojeve sa 100:

Take[data1,20]

data1=0.01{rmrf,foodrf}\[Transpose];

model1=LinearModelFit[data1,xrmrf,xrmrf]

Normal[model1]

Evo kako fit izgleda grafički

Show[ListPlot[model1["Data"],PlotStyle->Orange],Plot[model1[x],{x,-15,15}]]

Evo i ocenjenog modela kao i osnovnih njegovih karakteristika

Normal[model1]

model1["ParameterConfidenceIntervalTable"]

model1["AdjustedRSquared"]

Regresiona konstanta je alfa=0.0032, dok je beta=0.75. Beta određuje u kojoj meri promene uslova na tržištu utiču na promenu prinosa na vrednosno otežinjeni portfolio za hranu.

Dodavanje drugog faktora u regresiju

Ukoliko kretanje tržišta u celini hvala ceo rizik, to jest ukoliko bi CAPM važio u potpunosti, alfa bi morala biti jednaka nuli, što očigledno nije slučaj. Stoga zaključujemo da bismo za konstrukciju pravilnog modela morali da uključimo i neke druge faktore.

Jedan faktor koji može da ima uticaj jeste takozvani smb faktor, koji predstavlja dodatni prinos malih kompanija u odnosu na prinose velikih kompanija. Naime, poznato je da male kompanije u srednjem ostvaruju veće prinose od velikih, između ostalog zato što su rizičnije i, naravno, mogu lakše da rastu po velikim stopama.

Razmatrimo (opet moramo da podelimo sve brojeve sa 100) i ovaj faktor zajedno sa rmrf:

data2=0.01{smb,rmrf,foodrf}\[Transpose];

Take[data2,20]

model2=LinearModelFit[data2,{xsmb,xrmrf},{xsmb,xrmrf}]

Normal[model2]

model2["ParameterConfidenceIntervalTable"]

model2["AdjustedRSquared"]

Primetimo da dodavanjem samo jedne varijable smb jesmo malo povećali Adjusted (R^2) ali nismo smanjili konstantu. Sve variable su pojedinačno statistički značajne.

Dodajmo i treći faktor

Dodajmo sada i hml u regresionu jednačinu. To je prinos na portfolio koji je konstruisan da uzme u obzir razliku u prinosima između value i growth akcija, odnosno akcija koje imaju sporiji i brži rast. Sada regresija postaje:

data3=0.01{smb,hml,rmrf,foodrf}\[Transpose];

model3=LinearModelFit[data3,{xsmb,xhml,xrmrf},{xsmb,xhml,xrmrf}]

Normal[model3]

model3["ParameterConfidenceIntervalTable"]

model4["AdjustedRSquared"]

model3["AdjustedRSquared"]

Primetimo da smo opet dobili da su na 95 nivou pouzdanosti sve promenljive statistički značajne. Takođe, alfa je i dalje značajna.

Objašnjavajuće varijable mogu da ulaze u regresiju i nelinearno. Recimo, ukoliko želimo da pored linearnog razmatramo i kvadratni uticaj promenljive smb to možemo da uradimo ovako:

model4=LinearModelFit[data3,{xsmb,(xsmb^2),xhml,xrmrf},{xsmb,xhml,xrmrf}]

Normal[model4]

Vremenske serije

Uvod

Veliki broj funkcija je definisan specifično za analizu vremenskih serija. Prvo, treba da znamo kako da napravimo neku vremensku seriju. Jedan način jeste da unesemo jedan skup tačaka od kojih su prve koordinate vreme odnosno Times, a druge su vrednost odnosno Values. Komanda kojom definišemo da je nešto vremenska serija je TimeSeries. Evo primera:

prvaVremenskaSerija=TimeSeries[{{1,1},{2,4},{3,0},{4,3},{5,2}}]

DateListPlot[prvaVremenskaSerija]

TimeSeries[FinancialData["TSLA",{{2015}}]]//DateListPlot

prvaVremenskaSerija["Values"]

prvaVremenskaSerija["Times"]

Prisetimo se da smo vremenske serije ranije dobijali kako iz Monte Carlo simulacija stohastičkih procesa tako i koristeći TimeSeries na FinancialData podatke kada smo pravili portfolije.

Razmotrimo sada vremensku seriju dnevnih vrednosti SP500 indeks. Kao što smo rekli, SP500 indeks jeste indeks koji prati kretanje 500 najvećih kompanija u SAD. Njegova dinamika aproksimira kretanje američkog tržišta akcija.

SP500=TimeSeries[FinancialData["SP500",{{2000,1,1},{2010,1,1}}]]

Ovako rasklapamo i sklapamo vremensku seriju ukoliko želimo da je snimimo na Hard Disk i posle ponovo napravimo!

Ukoliko internet radi, ovaj deo možemo da preskočimo. Šta da radimo ako internet ne radi (ili na ispitu)? Jedan način da snimimo neku vremensku seriju jeste da je rekonstruišemo na “Times” i “Values” komponente. Onda svaku od njih možemo sačuvati kao .dat file. Pri tome, potrebno će nam biti na konvertujemo datume.

Snimili smo posebno "Times" i posebno "Values" na fajlove SP500Times2000do2010.dat i SP500Values2000do2010.dat respektivno pomoću Export funkcije. Sada možemo da iskoristimo te fajlove da rekonstruišemo originalnu SP500 vremensku seriju. Ovako ih snimamo na hard disk.

Export["SP500Times2000do2010.dat",SP500["Times"]]

Export["SP500Values2000do2010.dat",SP500["Values"]]

Primetimo da je potrebno da transformišemo vreme iz apsolutnog u kalendarsko vreme. Za prvu varijablu potrebni su nam podaci o vremenu, to jest file SP500Times2000do2010.dat

apsolutnoVreme=Import["C:\\Users\\Branko Urosevic\\Desktop\\OI 2016 Predavanja\\Podaci za WL Primene\\SP500Times2000do2010.dat"]\[Transpose][[1]];

datumi=DateObject/@apsolutnoVreme;

vrednosti=Import["C:\\Users\\Branko Urosevic\\Desktop\\OI 2016 Predavanja\\Podaci za WL Primene\\SP500Values2000do2010.dat"]\[Transpose][[1]];

SP500=TimeSeries[{datumi,vrednosti}\[Transpose]]

Analizirajmo sada SP500 vremensku seriju

DateListPlot[SP500]

Primetimo da je ovaj indeks bio na najnižem nivou u 2002, posle kraha dot.com balona, kao i u 2008-2009, u vreme velike svetske ekonomske krize. Naš naredni zadatak jeste da izračunamo logaritamske dnevne prinose na ovaj indeks. Primetimo da sve numeričke Listable funkcije možemo da koristimo direktno na vremenskim serijama. Konkretno, primenićemo prvo funkciju Log pa onda Differences. Evo kako to radi na nekoj apstraktoj seriji podataka:

Differences[Log[{x1,x2,x3,x4,x5}]]

Izračunajmo sada vremensku seriju logaritamskih prinosa za SP500. Kao rezultat opet dobijamo vremensku seriju, ali ovog puta logaritamskih dnevnih prinosa na akcije.

logprinosi=Differences[Log[SP500]]

DateListPlot[logprinosi,PlotRange-> All]

Dnevni prinosi osciluju oko nule. Takođe, dani velikih promena odnosno dani velike volatilnosti su grupisani a isto tako i dani kada je na tržištu sve mirno. Najveća volatilnost je bila, kao što vidimo, u vreme velike finansijske krize 2008-2009, kao i u periodu kada je propao dot.com balon.

Nije teško ustanoviti da se očekivana vrednost prinosa na dnevnom nivou vrlo malo razlikuje od nule. Naime:

Mean[logprinosi["Values"]]

LocationTest[logprinosi["Values"],0]

Modeliranje varijanse prinosa na akcije pomoću pokretnog proseka

Pošto je dnevni očekivani prinos praktično nula, ukoliko želimo da modeliramo kvadrat volatilnosti, odnosno varijansu možemo je bez problema aproksimirati pomoću kvadrata prinosa.

kvadprinosa=(logprinosi^2)

DateListPlot[kvadprinosa,PlotRange-> All]

U modeliranju vremenskih serija, u cilju boljeg sagledavanja trendova obično moramo da na neki način uprosečimo podatke. Jedan način jeste da napravimo pokretni prosek. Ukoliko imamo neku grupu podataka, pokretni prosek sa prozorom od 3 dana dobili bismo ovako:

MovingAverage[{x1,x2,x3,x4,x5},3]

Primenimo sada 100-dnevni pokretni prosek na kvadrat logaritamskih prinosa

MovingAverage[kvadprinosa,50]

Crvena linija predstavlja pokretni prosek od 100 dana

Show[DateListPlot[kvadprinosa,PlotRange-> All],DateListPlot[MovingAverage[kvadprinosa,100],PlotRange-> All,PlotStyle-> Red]]

Vidimo da 100-dnevni prosek u značajnoj meri poravnjava vremensku seriju.

Modeliranje varijanse prinosa na akcije pomoću EWMA

Ono što je problematično u vezi sa običnim pokretnim prosekom jeste da svi podaci u okviru prozora imaju isti ponder. U stvarnosti, prethodnih nekoliko dana trgovanja mnogo više utiče na tržište nego što je slučaj sa onim što se desilo pre 3 meseca. Stoga bi svežiji podaci morali da imaju veći ponder nego stari podaci. Jedan od najjednostavnijih modela koji omogućava da damo veći ponder svežijim podacima jeste takozvani ExponentialMovingAverage ili EWMA.

ExponentialMovingAverage[{x1,x2,x3,x4,x5},λ]

U praksi najčešće koristimo pondere od λ=0.06. To je, naime, parametar koji je uveo JP Morgan i konzorcijum velikih banaka kada su u opticaj uveli modele za pravljenje rizika, takozvane Value-at-Risk modele. Koristeći taj parametar predstavimo originalnu seriju podataka kao i varijansu modeliranu pomoću EWMA modela.

Show[DateListPlot[kvadprinosa,PlotRange-> All],DateListPlot[ExponentialMovingAverage[kvadprinosa,0.06],PlotRange-> All,PlotStyle-> Red]]

Vidimo da EWMA više prati stvarne promene volatilnosti bolje nego običan pokretni prosek od 100 dana. Takođe, za razliku od običnog pokretnog proseka koji kao rezultat daje vremensku seriju sa onoliko manje tačaka kolika je dužina perioda za koji se izračunava prosek (u našem slučaju 100 dana), kod EWMA to nije slučaj pa su originalna i EWMA serija iste dužine. EWMA model je jedan od najčešće korišćenik modela volatilnosti u modelima za upravljanje rizicima. Koristi ga veliki broj banaka uključujući, na primer, i Banca Intesa u Beogradu.

Modeliranje volatilnosti prinosa pomoću pokretnih proseka i EWMA

Primetimo da ukoliko nas zanima dnevna volatilnost a ne varijansa, pomoću pokretnog proseka i EWMA procenili bismo je prosto uzimanjem kvadratnog korena od vremenske serije za varijansu:

volatilnostprosek=(Sqrt[MovingAverage[kvadprinosa,100]])

volatilnostEWMA=√ExponentialMovingAverage[kvadprinosa,0.06]

Show[DateListPlot[volatilnostprosek,PlotStyle->Blue],DateListPlot[volatilnostEWMA,PlotStyle-> Red,PlotRange-> Full],PlotRange-> All]

Modeliranje vremenske serije za devizni kurs EUR-USD

Razmotrimo kako se na dnevnom nivou kretao kurs evra prema dolaru u 2015. Ovo poziva Wolfram Alpha:

kursEURpremaUSD=TimeSeries[First[WolframAlpha["euro to dollar in 2015",{{"Result",1},"ComputableData"}]]]["Values"];

Ako internet radi, onda ovo možemo da preskočimo. Međutim, ako ne radi kako treba, treba da snimimo vrednosti vremenske serije.

Pomoću ove komande ćemo očistiti od jedinica ovu vremensku seriju. Kao rezultat dobićemo samo brojeve, što onda lako možemo da sačuvamo.

kursEURpremaUSD=QuantityMagnitude[#]&/@kursEURpremaUSD;

Ovako možemo da eksportujemo ove vrednosti u fajl kursEURpremaUSD.dat

Export["kursEURpremaUSD.dat",kursEURpremaUSD]

A ovako možemo ponovo da učitamo vrednosti kursa u naš Notebook (ova linija koda nam treba ako nemamo pristup internetu):

kursEURpremaUSD=Import["C:\\Users\\mathematica\\Desktop\\Podaci za WL Primene\\kursEURpremaUSD.dat"]\[Transpose][[1]];

Sada smo spremni da radimo sa ovom veličinom.

Testiranje vremenske serije na jedinični koren

Pretpostavimo da nam je cilj da modeliramo devizni kurs EUR/USD koristeći AR, MA i ARMA modele. Devizni kurs dobijamo importujući fajl kursEURpremaUSD.dat

kursEURpremaUSD=Import["C:\\Users\\Branko Urosevic\\Desktop\\OI 2016 Predavanja\\Podaci za WL Primene\\kursEURpremaUSD.dat"]\[Transpose][[1]];

Evo kako izgleda vremenska zavisnost deviznog kursa sa vremenom:

ListLinePlot[deltadevizniKurs]

ListLinePlot[kursEURpremaUSD]

Proverimo da li ova vremenska serija ima jedinični koren. Probajmo, recimo:

jedinicnikorennivo=UnitRootTest[kursEURpremaUSD,1,"HypothesisTestData",SignificanceLevel->.001];

jedinicnikorennivo["TestConclusion","DickeyFullerT"]

Vidimo, stoga, da ova vremenska serija ima jedinični koren. Kao takva, ona nije pogodna za modeliranje pomoću ARMA modela (morali bismo da fitujemo ARIMA model u tom slučaju). Ukoliko prodiferenciramo jednom ovu seriju dobijamo:

deltadevizniKurs=Differences[kursEURpremaUSD];

jedinicnikorendelta=UnitRootTest[deltadevizniKurs,1,"HypothesisTestData",SignificanceLevel->.1];

jedinicnikorendelta["TestConclusion","DickeyFullerT"]

Modeliranje deviznog kursa korišćenjem AR, MA i ARMA procesa

Pretpostavimo da želimo da modeliramo devizni kurs kao proces AR(3) proces. AR proces dobijamo kada imamo zavisnost procesa od 1., 2., 3. i td lagova samog procesa, dok je greška izračunata u istom trenutku kao i proces sa leve strane. Drugim rečima, želimo da obasnimo proces na osnovu prethodnih realizacija procesa. Ukoliko je proces AR(3), koji ćemo mi ovde testirati, imali bismo:

(Subscript[fx, t])=a (Subscript[fx, t-1])+b (Subscript[fx, t-2])+c (Subscript[fx, t-3])+(Subscript[ϵ, t])

eAR3=EstimatedProcess[deltadevizniKurs, ARProcess[{a,b,c},v]]

MA proces dobijamo kada umesto lagova za vremensu seriju koju proučavamo koristimo lagove za grešku. U slučaju MA(3) koristimo 3 laga za grešku.

(Subscript[fx, t])=(Subscript[ϵ, t])+a (Subscript[ϵ, t-1])+b (Subscript[ϵ, t-2])+c (Subscript[ϵ, t-3])

Najzad, ARMA modeli kombinuju AR i MA procese. Mi ćemo koristiti ARMA(2,1) model gde se diferenca deviznog kursa regresira sa dva laga, a greška sa jednim. Ovi procesi su:

(Subscript[fx, t])=a1 (Subscript[fx, t-1])+a2 (Subscript[fx, t-2])+(Subscript[ϵ, t])+b (Subscript[ϵ, t-1])

Ove procese možemo da procenimo na sledeći način:

eMA3=EstimatedProcess[deltadevizniKurs, MAProcess[{a,b,c},v]]

eARMA21=EstimatedProcess[deltadevizniKurs, ARMAProcess[{a1,a2},{b},v]]

Upoređenje koeficijenata autokorelacije

Sada možemo da uporedimo prvih 10 koeficijenata autokorelacije za originalni devizni kursa kao i za procenjene procese:

ρdeviznikurs=Rest[CorrelationFunction[deltadevizniKurs,{11}]];

ρAR3=Rest[CorrelationFunction[eAR3,{11}]["Values"]];

ρMA3=Rest[CorrelationFunction[eMA3,{11}]["Values"]];

ρARMA21=Rest[CorrelationFunction[eARMA21,{11}]["Values"]];

Najzad, možemo da istom grafiku da predstavimo sve korelacije:

ListLinePlot[{ρdeviznikurs,ρAR3,ρMA3,ρARMA21},DataRange->{0,10},PlotLegends->{"EURUSDkurs","AR3","MA3","ARMA21"},PlotRange-> All]

Primetimo da kod promena deviznog kursa ono što se desilo sa promenom kursa pre deset data definitivno i dalje utične na promenu kursa danas. Vizualno liči da ARMA(2,1) najbolje odgovara podacima, mada daleko od savršenog. Najmanje odgovara MA(3).

Slaba stacionarnost stohastičkih procesa

Jedna od važnih karakteristika stochastičkih procesa jeste da li su slabo stacionarni. To znači prva dva momenta nisu funkcije vremena. Primetimo da su sva tri modela stacionarna u ovom smislu.

WeakStationarity[#]&/@{eAR3,eMA3,eARMA21}

Sa druge strane, na primer, WienerProcess nije stacionaran u ovom smislu zato što varijansa se menja kao funkcija vremene (linearno).

WeakStationarity[WienerProcess[0,σ]]

WeakStationarity[WienerProcess[μ,σ]]

Nije teško ustanovi ni da je PoissonProcess nestacionaran. Naime kod njega i srednja vrednost i varijansa se menjaju (linearno) sa vremenom:

Mean[PoissonProcess[λ][t]]

Variance[PoissonProcess[λ][t]]

Stoga je

WeakStationarity[PoissonProcess[λ]]