Logistinen regressio ja luottoriskien estimointi R-ohjelmalla

Tässä tekstissä käsittelen logistista regressiota ja luottoriskien estimointia. Muodostan mallin, jolla pankin asiakkaiden maksukyvyttömyyden (default) todennäköisyyttä voidaan arvioida perustuen erilaisiin muuttujiin. Muuttujina ovat mm. asiakkaan luottohistoria, haetun luoton määrä, työllisyystilanne, säästöt, asumismuoto ja ikä. Data (german.data) ja lisätietoa datasta löytyy täältä. Menetelmät, joita käytän perustuvat D, Sharman oppaaseen, joka löytyy täältä. Datasta Default-muuttuja toimii selitettävänä muuttujana ja muut 20 muuttujaa selittävinä muuttujina.
Default-muuttuja on binäärimuuttuja, joka saa arvon 1 jos asiakas on maksukyvytön ja arvon 0 jos asiakas ei ole maksukyvytön, eli pystyi maksamaan luottonsa takaisin. Tarvitsemme siis luotonmyöntöön vaikuttavien tietojen lisäksi tiedon siitä, kuinka asiakkaille kävi, eli pystyivätkö he maksamaan luottonsa vai eivät. Tämän tiedon avulla voimme muodostaa mallin ja mahdollisesti tehdä havaintoja, kuten esim. että työssäkäyvät selviytyvät luotoistaan työttömiä useammin. Jos nyt kaksi muuten täysin identtistä asiakasta tulisi hakemaan pankista lainaa, niin mallimme mukaan työllinen selviytyisi luotostaan todennäköisemmin. Se, että toinen maksaa lainansa takaisin todennäköisemmin ei tarkoita välttämättä, että toinen ei saisi lainaa ollenkaan. Lainanantaja voi nimittäin korkoa kasvattamalla pitää tuotto-odotuksensa vakiona, eli vaikka maksukyvyttömyyttä esiintyisi enemmän tietyllä ryhmällä, niin suuremmat korot voivat kompensoida nämä tappiot. Alla on esimerkki yksinkertaisesta yhtälöstä, joka kuvaa riskillisen tuoton ja riskittömän tuoton yhteyttä. Riskitön tuotto "r(g)" voi kuvata esim. EKP:n talletuskorkoa, eli tuottoa, jonka pankki saa rahalleen ilman riskiä. Yhtälössä "(1-p)" on "ei maksukyvyttömyyden" todennäköisyys, "u" on raha, joka saadaan takaisin maksukyvyttömyyden sattuessa (oletetaan, että kaikki menetetään, jos asiakas on maksukyvytön, eli u=0), "z" riskipreemio ja "r(f)" riskillinen tuotto . Riskipreemio kuvaa tuottoa, jonka pankki haluaa lisäksi, kun riskillisen tuoton odotusarvo vastaa riskitöntä tuottoa. Yhtälöstä (1) voidaan ratkaista ns. spreadi, joka kuvaa riskillisen- ja riskittömän tuoton erotusta (S=r(f)-r(g)). Tämä spreadi on esitetty yhtälössä (2). (Huom. ratkaisussa on oletettu, että p ja r(f) ovat pieniä, jotta p*r(f)~0

1.  (1-p)(1+r(f)) + p*u = 1+r(g)+z
2.  S=p+z

Spreadiin vaikuttaa siis sekä maksukyvyttömyyden todennäköisyys että riskipreemio. Vielä yksinkertaistan, kun oletetaan, että riskipreemio on nolla, voimme kirjoittaa spreadin S=p eli r(f)-r(g)=p

Jos siis pankki ei pyydä riskipreemiota ja haluaa odotusarvoisesti saada saman tuoton sekä työlliseltä että työttömältä, niin sen pitää tietää p mahdollisimman tarkasti. Jos työllisen todennäköisyys maksukyvyttömyyteen olisi 1%:n ja työttömän vastaava luku 3%:a, niin pankki pyytäisi esimerkin tilanteessa työlliseltä korkoa ¨~1%+riskitön korko ja vastaavasti työttömältä  ~3%+riskitön korko. Tällä tavoin pankin odotusarvo kummallekin saatavalla olisi sama.


Nyt takaisin malliin ja sen luomiseen. Ensin lataan R:ään datan ja nimeän sarakkeet, koska niillä ei ole valmiina nimiä.
R-koodi:

rm(lists=ls())    #poistaa aiemmat tallennetut muuttujat. Hyvä käyttää, jos attach-komentoa käytetään.
data <- read.delim(file.choose(), header=FALSE,sep="")

names(data) <- c("check_acc", "duration", "history", "purpose", "amount", "sav_acc", "employment", "install_rate", "ptatus", "other_debtor", "resid", "property", "age", "other_install", "housing", "other_credits", "job", "num_depend", "teleph", "foreign", "default")
attach(data)
str(data)                              #str-komennolla nähdään ovatko muuttujat kategoriallisia vai jatkuvia

Seuraavaksi voimme katsoa dataa kuvien avulla. Jatkuvia muuttujia, kuten age voidaan katsoa esim. plot-ja hist-komennoilla. Muuttujia, jotka ovat faktoreina (esim.history) voidaan katsoa lattice-pakkauksen histogram-komennolla. Katsotaan seuraavaksi age-muuttujaa:
R-koodi:

library(lattice)

histogram(age,type="percent",breaks=seq(10,80,by=10))   #breaks-komennolla jaetaan eri ikäryhmiin

Näemme, että aineistossa suurin osa on 20-40-vuotiaita. Voimme tarkastella eri ikäluokkien osuuksia default-luokissa vcd-pakkauksen mosaic-komennolla. Jaetaan age-muuttuja ensin kolmeen ikäluokkaan ja muutetaan default-muuttuja faktoriksi, joka saa arvot 0 ja 1. Sen arvo 1 tarkoittaa maksukyvyttömyyttä. Sitten voimmekin tehdä Mosaic-kuvan.
R-koodi:

agec <- as.factor(ifelse(age <= 30, "0-30",ifelse(age <=40,"31-40","41-80")))
data$default <- as.factor(ifelse(data$default==1,"0","1")) #päivitetään myös data.frameen
attach(data)
library(vcd)
mosaic(~default+age,highlighting="age",main="age & default")

Kuvasta nähdään, että 0-30-vuotiaiden osuus on maksukyvyttömistä suurin. Seuraavaksi tutkin, mitkä ovat maksukyvyttömien osuudet eri ikäluokissa. Tämän voi tehdä table-komennolla seuraavasti: 
R-koodi:

defaultagetable <- table(default,age) # Defaultin luokat menevät table:n riveiksi ja age:n luokat sarakkeiksi
prop.table(defaultagetable,margin=2) # Lasketaan maksukyvyttömien ja maksukykyisten osuudet
a <- barplot(prop.table(defaultagetable,margin=2),beside=T,legend.text=c("Maksukykyisten osuus",
"Maksukyvyttömien osuus"),yaxt="n",ylab="%",xlab="age") # Tehdään pylväskuvio
axis(2,at=seq(0,1,by=0.1))
text(a, 0, round(fdf,2),cex=1.1,pos=3)

 Testaan nyt Khiin neliö -testillä (Chi-square test), ovatko nämä erot ikäluokkien välillä tilastollisesti merkitseviä. "0-30" -luokka eroaa muistaa luokista merkitsevästi, kun taas "31-40" ja "41-80" -luokat eivät eroa keskenään merkitsevästi.

Seuraavaksi lisään "age & default" -kuvaan kolmannen muuttujan, joka on other_debtor, eli muuttuja, joka on kuvattu aineistokuvauksessa seuraavasti:

Other debtors / guarantors 
A101 : none 
A102 : co-applicant 
A103 : guarantor 


Muuttuja siis kertoo, onko hakijalla lainalle takaaja ja hakeeko hakija lainaa yksin. Katsotaan seuraavaksi kuvaa.
R-koodi:

mosaic(default~age+other_debtor,main="age, default & other_debtor")

Kuvasta nähdään, että jokaisessa ikäryhmässä "co-applicant"-hakijat olivat suhteellisesti useimmin maksukyvyttömiä. Hakijat, joilla oli takaaja, pystyivät joka ikäluokassa suhteellisesti useimmin maksamaan lainansa takaisin. 

Seuraavaksi muutan osan muuttujista faktori-muotoon, koska ne eivät ole numeerisia muuttujia eli ne saavat arvoja esim. välillä 1-4 ja kuvaavat luokittelua. 
R-koodi:

data$install_rate <- as.factor(data$install_rate)
data$resid <- as.factor(data$resid)
data$other_credits <- as.factor(data$other_credits)

data$num_depend <- as.factor(data$num_depend)

Nyt  jaan datan "trainingdataan" ja "testdataan", eli jaan datan kahteen osaan, joista toisesta (trainingdata) muodostan mallin, jolla ennustan toisen arvoja (testdata). Tosiasiassa mallin toimivuus tulisi testata myös cross-validationilla ja/tai jakamalla data kolmee osaan, joista kolmatta validationsettiä ei käytetä mallia luodessa, vaan sitä kokeillaan vain kerran, kun malli on valmis. Kun mallia verrataan testdataan, niin muuttujia lisäillessä ja mallia "tuunatessa" voi syntyä ns. overfitting, jolloin malli toimii testdataan paremmin kuin mallin ulkopuoliseen validationsettiin. Tavoitteena on siis muodostaa malli, johon syöttämällä testdatan muuttujien arvot voidaan saada mahdollisimman hyviä ennusteita maksukykyisyydestä. Testdatastahan voimme tarkistaa, kuinka hyvin malli on onnistunut, koska testdatassa on default-muuttuja, joka kertoo maksoiko asiakas lainansa takaisin vai ei (tosiasiassa 1-arvo voi tarkoittaa myös, että asiakas on maksanut lainansa myöhässä tai lainanmaksussa on ollut muita vaikeuksia). Otan trainingdataan sattumanvaraisesti 60%:a arvoista, eli testdataan jää 40%:a arvoista. Huom. Datassa on tällä hetkellä mukana myös jatkuvia muuttujia, kuten age ja amount, joten teen aluksi mallin normaalilla datalla. Kokeilen tämän jälkeen mallia, jossa on pelkästään kategoriallisia muuttujia, eli myös muuttujat, kuten age ja amount on muutettu kategorioihin. Sharma suosittelee käytettäväksi mallissa vain kategorisia muuttujia, mutta toisaalta informaatiokadon vuoksi myös jatkuvien muuttujien malleilla on kannattajansa. Jatkuvan muuttujan muuttaminen kategorialliseksi voi olla hyödyllistä, jos muuttujan yhteys selitettävään muuttujaan ei ole lineaarinen. Epälineaarisia yhteyksiä voi toisaalta mallintaa myös esim. neliöjuuri-muunnoksella, jolloin mallissa käytetään muuttujaa sekä normaalissa että neliöjuuri-muodossa. Epälineaarisuuksia tarkastellessa on aina tarkasteltava muuttujaa myös eri kategorioissa, koska epälineaarisuus saattaa johtua  muuttujan erilaisesta käyttäytymisestä jossain kategoriassa (yhteyden ollessa lineaarinen muissa kategorioissa). Jos muuttuja käyttäytyy eri tavoin eri kategorioissa, on interarction-termien käyttö suositeltavaa. Suosittelen, että muuttujia tarkastellaan aina ennen mallin tekemistä huolella ja mahdollinen kategorisointi on pystyttävä perustelemaan teorialla. Lähtökohtaisesti pyrkisin pitämään jatkuvan muuttujan jatkuvana. Joskus selitettävän ja selittävän -muuttujan välillä löytyy yhteys vasta esim. logaritmi-muunnoksen jälkeen, joten eri muunnokset on hyvä käydä läpi ennen kuin lähtee kategorisoimaan muuttujaa.
R-koodi:

e <- sample(nrow(data), 0.6*nrow(data))  #ottaa sattumanvaraisesti 60%:a rivinumeroista
trainingdata <- data[e,]                              # valitaan saadut rivinumerot trainingdataan ja testdataan

testdata <- data[-e,]

Nyt trainingdatasta voidaan muodostaa logistinen regressiomalli. 
R-koodi: 

model1 <- glm(default~check_acc+duration+history+purpose+
amount+
sav_acc+
employment+
install_rate+
ptatus+
other_debtor+
resid+
property+
age+
other_install+
housing+
other_credits+
job+
num_depend+
teleph+

foreign,family=binomial(link=logit),data=trainingdata)
summary(model1)

Mallissa 1 on useampia tilastollisesti merkitseviä muuttujia. Alla tuloksia. (Standard error suluissa)

Results of the model1

	Dependent variable:
	default
check_accA12	-0.550*
	(0.282)
check_accA13	-0.933**
	(0.461)
check_accA14	-1.944***
	(0.321)
duration	0.036***
	(0.012)
historyA31	0.849
	(0.727)
historyA32	-0.059
	(0.565)
historyA33	-0.240
	(0.613)
historyA34	-0.527
	(0.557)
purposeA41	-1.859***
	(0.496)
purposeA410	-2.696**
	(1.134)
purposeA42	-0.553
	(0.346)
purposeA43	-0.809**
	(0.319)
purposeA44	-2.336
	(1.480)
purposeA45	-0.166
	(0.650)
purposeA46	0.092
	(0.563)
purposeA48	-0.696
	(1.465)
purposeA49	-0.297
	(0.423)
amount	0.0001
	(0.0001)
sav_accA62	-0.541
	(0.377)
sav_accA63	-0.532
	(0.502)
sav_accA64	-1.364*
	(0.734)
sav_accA65	-1.044***
	(0.365)
employmentA72	0.759
	(0.637)
employmentA73	0.454
	(0.617)
employmentA74	0.070
	(0.654)
employmentA75	0.517
	(0.614)
install_rate2	0.150
	(0.389)
install_rate3	0.332
	(0.439)
install_rate4	0.731*
	(0.387)
ptatusA92	-0.084
	(0.490)
ptatusA93	-0.630
	(0.477)
ptatusA94	0.296
	(0.579)
other_debtorA102	0.331
	(0.637)
other_debtorA103	-0.844*
	(0.505)
resid2	1.186***
	(0.400)
resid3	0.872**
	(0.429)
resid4	0.452
	(0.403)
propertyA122	0.446
	(0.325)
propertyA123	0.139
	(0.312)
propertyA124	0.950*
	(0.554)
age	-0.019
	(0.012)
other_installA142	-0.682
	(0.651)
other_installA143	-0.520
	(0.318)
housingA152	-0.391
	(0.314)
housingA153	-0.221
	(0.636)
other_credits2	0.455
	(0.312)
other_credits3	0.073
	(0.848)
other_credits4	0.851
	(1.384)
jobA172	-0.169
	(0.807)
jobA173	-0.198
	(0.779)
jobA174	0.257
	(0.807)
num_depend2	0.261
	(0.328)
telephA192	-0.372
	(0.265)
foreignA202	-2.272***
	(0.854)
Constant	-0.051
	(1.309)
Observations	600
Log Likelihood	-270.011
Akaike Inf. Crit.	650.022
Note:	*p<0.1; **p<0.05; ***p<0.01

>

Nyt voimme arvioida mallin ennustustarkkuutta laskemalla ROC-käyrän.
R-koodi:

testdata$model1p <- predict(model1, testdata, type = "response") # predict-komento laskee mallilla ennusteet testdatan arvoilla ja muodostaa testdataan uuden "model1p"-sarakkeen
library(ROCR)
model1points <- prediction(testdata$model1p, testdata$default) #vertaa ennusteita oikeisiin arvoihin
performancemodel1 <- performance(model1points,"tpr","fpr")
plot(performancemodel1, col = "blue")     #kuvaa ROC-käyrän
abline(0,1,lty=5)                                         #tekee 45 asteen suoran viivan

ROC-käyrän muotoja ja mallin onnistumistarkkuutta havainnollistaa alla oleva kuva.

Lähde: http://gim.unmc.edu/dxtests/roccomp.jpg

ROC-käyrässä True positive rate (TPR) tarkoittaa lukua, joka saadaan, kun jaetaan onnistuneet luotonhylkäykset, kun luotonsaaja on maksukyvytön, kaikilla maksukyvyttömien määrällä. False positive rate (FPR) tarkoittaa lukua, joka saadaan, kun jaetaan luoton hylkäykset, kun asiakas olisi ollut maksukykyinen, kaikkien maksukykyisten määrällä. Käyrän arvot saadaan, kun käytetään erilaisia ns. cut-off -arvoja ja lasketaan TPR ja FPR näillä arvoilla. Esim. kun cut-off-arvo on 0.1, niin luotot hylätään, jos asiakkaan maksukyvyttömyyden todennäköisyys on yli tämän. Hylätyistä luotoista lasketaan sitten oikein ja väärin luokitellut, jolloin saadaan ROC-käyrän piste tälle cut-off-arvolle.  ROC-käyrästä näemme, että cut-off arvo alkaa luvusta 1 ja oikealle mentäessä se laskee aina nollaan saakka, jolloin käyrät kohtaavat. Lasken vielä AUC-arvon, joka kertoo, kuinka paljon aluetta ROC-käyrän alle jää. Tämä auttaa ROC-käyrien vertailussa, kun ero ei visuaalisesti ole aivan selvä. ROC-käyrille saa laskettua myös luottamusvälit. 
R-koodi:

model1_auc <- performance(model1points, "auc")

Mallin AUC-arvo on 0.7898

Lasken myös Kolmogorov-Smirnov (KS) -arvon
R-koodi:


ks1 <- max(attr(performancemodel1, "y.values")[[1]] - (attr(performancemodel1 "x.values")[[1]]))



KS-arvo kertoo FPR:n ja TPR:n kumulatiivisen erotuksen maksimikohdan. KS-arvo 1 tarkoittaa, että malli on onnistunut jakamaan maksukyvyttömät ja maksukykyiset selkeästi kahteen eri luokkaan. KS-arvo 0 tarkoittaisi, ettei malli ole onnistunut luokkiin jakamisessa. KS-arvo 1 voitaisiin saada esim. jos kaikki maksukykyiset saisivat arvon <0.5 ja kaikki maksukyvyttömät arvon >= 0.5. Tällainen KS-arvo tarkoittaisi, että voisimme laittaa cut-off arvoksi 0.5, jolloin hylkäisimme kaikki maksukyvyttömien lainat ja myöntäisimme kaikki maksukykyisten lainat. 

KS-arvosta löytyy lisää tietoa mm. täältä. 

Kokeilin myös mallia, jossa jatkuvat muuttujat (age, duration ja amount) oli muutettu kategorisiksi, mutta mallin tulokset eivät eronneet merkittävästi. Mallia sai paranneltua, kun lisäsi siihen interaction-muuttujia. Potentiaalisia interaction-muuttujia löytää Machine learningiin (Koneoppiminen) kuuluvan Random Forestin Variable importance -listasta. Parhaan mallin löysin käyttämällä Koneoppimisen menetelmiä, joilla sain AUC-arvon, joka oli suurempi kuin D. Sharman paperissaan saama suurin arvo. Koneoppiminen on erittäin hyödyllinen menetelmä mallintamisessa, koska sen keinoin saadaan useista malleista luotua ennusteita, jotka peittoavat yksittäisten mallien tulokset. Näistä menetelmistä kirjoittelen jatkossa tarkemmin. 

Aikasarjan autokorrelaatio

Autokorrelaatiokuvat ACF ja PACF ovat olennaisessa osassa, kun tutkitaan aikasarjoja ja halutaan ennustaa. On siis tärkeää ymmärtää, miten kuvissa esiintyvät luvut syntyvät. Katsotaan nyt myyntien aikasarjaa, jota käytin SARIMA-esimerkissä ja tutkitaan autokorrelaatiota tarkemmin. Olen lisännyt kuvaan viivan, joka kuvaa myyntien keskiarvoa.

Alla myös aikasarjan autokorrelaatiokuvat.

Autokorrelaatio lasketaan alla olevalla kaavalla. 

Autokorrelaation suuruuteen vaikuttaa etenkin y(t):n ja y(t-1) välinen kovarianssi. Ensimmäistä kuvaa katsomalla näemme, että noin puolet arvoista on keskiarvon alapuolella ja noin puolet keskiarvon yläpuolella. Tästä johtuen, kun laskemme y(t):n autokorrelaatiota ensimmäiselle viiveelle, niin y(t):n ja y(t-1):n kovarianssi on jatkuvasti positiivinen. y(t):n varianssi on luonnollisesti positiivinen, joten tällöin myös autokorrelaatio on ensimmäiselle viiveelle positiivinen. 

Annan esimerkin autokorrelaation laskemisesta ensimmäiselle viiveelle. Ajatellaan, että meillä on aikasarja, joka muodostuu arvoista 1,2,3 ja 4, ja aikaväli on 01/2015-04/2015 (eli tammikuun arvo on 1, helmikuun arvo on 2, maaliskuun arvo on 3 ja huhtikuun arvo on 4). Olemme kiinnostuneita siitä korreloiko aikasarja ensimmäisen viiveensä kanssa, eli selittääkö tammikuun arvo helmikuun arvoa ja helmikuun arvo maaliskuun arvoa ja maaliskuun arvo huhtikuun arvoa. Lasketaan ensin kovarianssi y(t:):lle ja y(t-1):lle. Aikasarjan keskiarvo on (1+2+3+4)/4= 2.5 joten voimme kirjottaa kovarianssin seuraavasti: (2-2.5)(1-2.5)+(3-2.5)(2-2.5)+(4-2.5)(3-2.5)= 1.25

Lasketaan myös varianssi: (1-2.5)^2 + (2-2.5)^2 + (3-2.5)^2 + (4-2.5)^2 = 5

Autokorrelaatio ensimmäiselle viiveelle on siis 1.25/5=0.25

Lasketaan autokorrelaatiot myös viiveille 2 ja 3. Ensimmäiseksi lasketaan kovarianssi 2. viiveelle (huom.varianssi pysyy samana viiveistä riippumatta): (3-2.5)(1-2.5)+(4-2.5)(2-2.5)= -1.5 eli autokorrelaatio 2. viiveelle on: -1.5/5= -0.3

Ja lopuksi 3. viiveen kovarianssi: (4-2.5)(1-2.5)= -2.25 eli autokorrelaatio 3. viiveelle on: -2.25/5= -0.45

Positiivinen autokorrelaatio tarkoittaa sitä, että havainnot valitulla viivellä ovat suureksi osaksi saman merkkisiä, eli joko keskiarvoa pienempiä tai suurempia. Autokorrelaatiota aiheuttavat erityisesti mallista puuttuvat muuttujat sekä vääränlaisen mallin käyttö. Puuttuva muuttuja voi olla esim. AR(1) termi, jota emme ole vielä lisänneet malliin. Vääranlaisella mallillla tarkoitan tässä mallia, jossa on vielä yksikköjuuri. Kuten SARIMA-esimerkissä näimme, niin aikasarjassa on ennen kausidifferensointia vahva autokorrelaatio, joka vähenee huomattavasti tämän toimenpiteen jälkeen. Toisaalta kaikki autokorrelaatio ei katoa mallista ja tarvitsemme puuttuvat muuttujat malliin. On tärkeää tiedostaa, että ARIMA-ennustaminen toimii, koska aikasarjoissa eri havainnot eivät usein ole täysin itsenäisiä. Eli aiemman kuukauden havainto tarjoaa usein informaatiota tulevan kuukauden havainnosta. Itseasiassa "weakly" stationaarisessa AR(1) mallissa ensimmäisen viiveen autokorrelaatio on sama kuin mallissa esiintyvän y(t-1) kerroin, eli ensimmäinen viive vaikuttaa y(t):n aina samalla tavalla. Jos havainnot olisivat täysin itsenäisiä ja riippumattomia toisistaan, eli ns.white noise, niin emme voisi ennustaa kuin, että aikasarja pysyy keskiarvossaan. (White noisella tarkoitetaan sarjaa, jonka havainnot pyörivät keskiarvonsa ympärillä täysin sattumanvaraisesti, joten aikaisemmat arvot eivät tarjoa informaatiota tulevasta). Tämä on hyvä pitää mielessä, jos haluaa ennustaa esim. osakemarkkinoita. Tarjoavatko osakkeen historialliset arvostukset oikeasti relevanttia tietoa tulevista kursseista vai onko paras ennustus vain viimeisin kurssi?

PACF

Osittaisautokorrelaatio mittaa lineaarista riippuvuutta tietyllä viiveellä, kun aikaisempien viiveiden vaikutus on otettu pois. Tästä johtuen aikasarjan autokorrelaatio ensimmäisellä viiveellä on sama kuin osittaisautokorrelaatio ensimmäisellä viiveellä. Osittaisautokorrelaatiokuvat kertovat kannattaako AR-termejä laittaa ARIMA-/SARIMA -malleihin. Osittaisautokorrelaatio siis kertoo, onko esim. aikasarjan 3.viiveellä (y(t-3) lineaarista riippuvuutta y(t):n kanssa. Jos näin ei ole, niin aikasarjamalliin ei kannata laittaa AR(3) termiä.

Osittaisautokorrelaatio kuvataan usein, niin että osittaisiautokorrelaatio viiveellä h olisi sama asia kuin kyseisen viiveen regressiokerroin mallissa, jossa ovat sen lisäksi sitä edeltävät viiveet. AR-malli olisi siis viiveen h=3 tapauksessa seuraava: y(t)=b0+b1y(t-1)+b2y(t-2)+b3y(t-3) ja osittaisautokorrelaatio viiveelle 3 olisi kerroin b3.  Vastaavasti 4.viiveen AR-malli olisi y(t)=b0+b1y(t-1)+b2y(t-2)+b3y(t-3)+b4y(t-4) ja sen osittaisautokorrelaatio olisi kerroin b4. Tätä menetelmää ei kuitenkaan R:ssä käytetä, vaan kyseessä on Durbin-Levinson Rekursio, joka on algoritmi, jolla vastaavat kertoimet (osittaisautokorrelaatiot) lasketaan. Tästä laskentatavasta voit lukea lisää esim. täältä

Seuraavaksi simuloin R:llä ei-stationaarisen aikasarjan, jotta näemme, kuinka perinteisen AR-mallin kertoimet ja R:n osittaisautokorrelaatioluvut eroavat toisistaan huomattavasti. Simuloin myöhemmin myös stationaarisen aikasarjan ja näemme, että AR-mallin kertoimet ovat hyvin lähellä osittaisiautokorrelaatiolukuja.

R-koodi ei- stationaariselle aikasarjalle:

set.seed(1)                                                    #set.seed(1):lla saat samat tulokset kuin tässä
x <- w <- rnorm(500)
for (t in 2:500) x[t]<-x[t-1]+w[t]                 # Tehdään random walk-aikasarja
par(mfrow=c(2,1))                                       #saadaan kaksi kuvaa samaan ruutuun
plot(x,type="l")
Pacf(x,lag.max=3)

Pacf(x,plot=FALSE,lag.max=3)                #saadaan tulostettua arvot kolmelle viiveelle.

Kuva (ylempi) ei-stationaarinen aikasarja. Kuva (alempi) Osittaisautokorrelaatio

Partial autocorrelations of series ‘x’, by lag

     1        2          3
 0.972     0.009   -0.023

Yllä osittaisautokorrelaatioluvut. Nyt lasken R:n Arima-funktiota käyttäen vastaaville viiveille arviot, eli teen AR-mallit 1.viiveelle, 2.viiveelle ja 3.viiveelle.
R-koodi:

Arima(x,order=c(1,0,0))                                       #Saadaan 1.viiveen osittaisautokorrelaatio
Arima(x,order=c(2,0,0))                                       #Saadaan 2.viiveen osittaisautokorrelaatio
Arima(x,order=c(3,0,0))                                       #Saadaan 3. viiveen osittaisautokorrelaatio

Näin saadut osittaisautokorrelaatiot ovat

1          2           3
0.979   0.013   -0.007

Näemme, etteivät arviot ole kovin lähellä toisiaan.

Nyt muodostan stationaarisen aikasarjan ja teen samat toiminnot kuin ei-stationaariselle, eli kuva aikasarjasta ja osittaisautokorrelaatiosta sekä osittaisautokorrelaatioluvut.
R-koodi:

set.seed(1)
w<-rnorm(500)
par(mfrow=c(2,1))
plot(w,type="l")
Pacf(w,lag.max=3)
Pacf(w,plot=FALSE,lag.max=3)

Kuva (ylempi) stationaarinen aikasarja. Kuva (alempi) Osittaisautokorrelaatio

Saadut osittaisautokorrelaatioluvut ovat:

     1         2           3
-0.027    -0.004   -0.040                                                              

Nyt muodostan Arima-funktiota käyttäen AR-mallit, kuten aiemmin. 

R-koodi:

Arima(x,order=c(1,0,0))

Arima(x,order=c(2,0,0))

Arima(x,order=c(3,0,0))

Nyt saadut osittaisautokorrelaatioluvut ovat:

1            2             3

-0.027  -0.004     -0,040

Luvut ovat siis pyöristettyinä samat. Stationaarisessa aikasarjassa osittaisautokorrelaatiot ovat siis lähelle AR-mallin ennustamia, mutta eivät aivan samoja. Esim. 5.viiveelle lasketut arvot ovat 0.022 (Pacf) ja 0.026 (AR(5)). Intuitiivisesti osittaisautokorrelaation voi kuitenkin ymmärtää juuri AR-mallin tavoin. Ajatellaan esimerkiksi kahden mallin eroa, kun haluamme laskea 2.viiveen autokorrelaation. 1. malli: y(t)=y(t-2) ja 2. malli y(t)=y(t-1)+y(t-2)

1.mallissa emme ota huomioon 1.viiveen vaikutusta, kun taas toisessa otamme. Koska haluamme tietää vain 2.viiveen vaikutuksen y(t):hen, pitää y(t-1) ottaa mukaan, jotta sen vaikutus saadaan myös huomioitua.
























Kirjoittaja on rekrytointisivusto: www.dataverkosto.com perustaja.