## MONIULOTTEISET AIKASARJAT ## R-KOODEJA VAR-MALLILLE; ks.lisää: http://www.jstatsoft.org/v27/i04/paper ## Tarvittavia paketteja install.packages("vars") ##var-mallin estimoinnissa library(vars) library(stats) ## Esimerkkinä käytetään aineistoa, jonka muuttujina ovat Länsi-Saksan kulutus, investoinnit ja tulot vuosilta 1960Q1-1982Q4 ## AINEISTON HAKU JA AINEISTON KÄSITTELY (voit olla tallentanut haluamasi tiedoston minne vain) ## header = true olettaa, että tiedoston ensimmäisellä rivillä on muuttujien nimet. Data = read.table(file.choose(), header=TRUE) ## Tämän komennon pitäisi antaa sinun valita kansioistasi haluamasi tiedoston. ## Vaihtoehtoinen tapa on. Muuta tiedoston polku oikeaksi ja poista kommentointi jos haluat käyttää tätä menetelmää. Data = read.table("C:/Opiskelu/Tilastotiede/Moniulotteiset 2015/Laskuharjoitukset/ger.txt",header=TRUE) attach(Data) ## Laittaa aineiston nimet tiedostokantaan (kumpaa tahansa edellä mainittua käytettäessä). ## Aineisto aikasarjaksi ## Esim. deltat = 1/12 olisi kuukausidataa. Seuraavassa neljännesvuosidataa. Frequency määritelmää voi myös käyttää. ger = ts(Data[,], start = c(1960,1), deltat = 1/4) ## Aikasarjojen piirtäminen plot.ts(ger) ## Tarvittaessa aineiston aikasarjojen logaritmointi ja differensointi (tässä molemmat). ld_ger = diff(log(ger)) plot.ts(ld_ger) ## Korrelaatiorakenteen tutkiminen ## (Yhtälöissä hieman manipulointia, jotta x-akseli tulisi järkeväksi. ## R sekoittaa jostain (tuntemattomasta) syystä x-akselin ja aikasarjan frekvenssin. ## Kun alempana asettaa ts(ld_ger,freq=1), niin tilanne korjaantuu. ## Tällä ei ole vaikutusta tuloksiin; kuvista tulee vain järkevämmät.) ## Auto- ja ristikorrelaatiofunktiot (vrt. monisteen Kuvio 3.1). ## Seuraavassa 'lag.max' määrittelee maksimiviipymän pituuden. acf(ts(ld_ger,freq = 1), lag.max = 10,type = "correlation") ## VAR-MALLIN (NORMAALIJAKAUMAAN PERUSTUVA EHDOLLINEN) SU-ESTIMOINTI ELI PNS-ESTIMOINTI ## Mallinvalintaa varten on koodi ?VARselect # Voit tarkastella koodin ominaisuuksia. ## Seuraavassa 'type' määrittelee estimoitavan mallin tyypin. ## Vakion (const) lisäksi voidaan sallia aikatrendi (trend), vakio ja aikatrendi (both) tai ei kumpaakaan (none) malval = VARselect(ld_ger, lag.max = 10, type = c("const"), season = NULL, exogen = NULL) ## Akaike (AIC), Hannan-Quinn (HQ), Schwarz (SC=BIC) ja Akaiken Final Prediction Error (viimeistä ei käsitellä kurssilla) malval$selection ## Tulostaa eri kriteerien ehdottamat asteet. ## VAR-mallin parametrien estimointi # p = mallin aste # type = onko mallissa vakio (const), aikatrendi (trend), molemmat (both), ei kumpaakaan (none) # Malliin voidaan sisällyttää myös muita selittäviä eli ns. eksogeenisia muuttujia (exogen) tai # kausivaihtelua huomioon ottavia indikaattorimuuttujia (eli 'kausidummeja'). # Jos eksogeenisia muuttujia käytetään, niistä on muodostettava matriisi. Tällaisia malleja ei käsitellä kurssilla, # joten tämä vain tiedoksi, samoin kuin se, että estimointitulosten yhteydessä analysoitavia muuttujia nimitetään endogeenisiksi muuttujiksi. # Seuraavassa asteeksi on valittu 2. var_lag2 = VAR(ld_ger, p = 2, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) # Voit tarkastella, mitä var_lag2 pitää sisällään str(var_lag2) summary(var_lag2) ## Estimointitulokset esitetään yhtälöittäin (ks. monisteen s. 22, yhtälöt (4.4)). # Estimoiduista yhtälöistä tulostetaan myös joitakin kurssilla käsittelemättömiä tietoja. # Otsikolla 'Roots of the characteristic polynomial' tulostetaan monisteen s. 8 esitetyn lihavoidun A-matriisin ominaisarvot. logLik(var_lag2) ## Laskee (ehdollisen) log-uskottavuusfunktion maksimiarvon. ## Mallinvalinta voidaan perustaa myös uskottavuusosamäärätesteihin, mutta valmista pakettia ei löytynyt, joten tämä on tehtävä manuaalisesti. ## Testi perustuu peräkkäiseen testaukseen. Oletetaan esimerkiksi, että VAR-mallin maksimaalinen aste on 4 ja lasketaan (ehdollisen) log-uskottavuusfunktion ## maksimiarvot (vaihtoehtoisesti testisuureiden arvot voitaisiin laskea virhekovarianssimatriisien estimaattien avulla kuten monisteessa s. 30). ## Jotta testaus menisi oikein, on varmistuttava siitä, että uskottavuusfunktion laskemisessa käytetään samaa määrää havaintoja (ks. monisteen s. 31). ## Tästä syystä eri VAR(p)-malleille on aineiston alkupäätä rajoitettava. Rakennetaan eri aineistoja siten, että jokainen estimoitava malli käyttää saman määrän havaintoja. ## Seuraavassa esim. start=c(1960,3) rajoittaa (logaritmoidun ja differensoidun sarjan ld_ger) havainnot alkamaan vuoden 1960 3. neljänneksestä. ## Tällöin estimoitavan VAR(3)-mallin selitettävän muuttujan ensimmäinen havainto on vuoden 1961 2. havainto, kuten maksimiviipymää p = 4 vastaavassa VAR(4)-mallissakin. ger_3 = window(ld_ger, start = c(1960,3)) ger_2 = window(ld_ger,start = c(1960,4)) ger_1 = window(ld_ger,start = c(1961,1)) ## VAR(4) var_lag4 = VAR(ld_ger, p = 4, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) logLik(var_lag4) ## VAR(3) var_lag3 = VAR(ger_3, p = 3, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) logLik(var_lag3) ## Uskottavuusosamäärätestisuureen arvo LR43 = 2*(logLik(var_lag4)-logLik(var_lag3)) ## Testissä käytettävä vapausasteiden lukumäärä määräytyy rajoitettavien parametrien lukumäärästä. ## Seuraavassa var_lag4$K on muuttujien lukumäärä (monisteen n) tässä sovellettavassa VAR(3)-mallissa eli kolme. ## Nyt testissä rajoitetaan nollaksi matriisi A_4, jossa on 3^2 parametria, mistä vapausasteiden lukumäärä määräytyy. p_arvo43 = 1-pchisq(LR43,(var_lag4$K)^2) p_arvo43 ## VAR(2) var_lag2 = VAR(ger_2, p = 2, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) logLik(var_lag2) ## Uskottavuusosamäärätestisuureen arvo LR32 = 2*(logLik(var_lag3)-logLik(var_lag2)) p_arvo32 = 1 - pchisq(LR32,(var_lag4$K)^2) p_arvo32 ## VAR(1) var_lag1 = VAR(ger_1, p = 1, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) logLik(var_lag1) ## Uskottavuusosamäärätestisuureen arvo LR21 = 2*(logLik(var_lag2)-logLik(var_lag1)) p_arvo21 = 1 - pchisq(LR21,(var_lag4$K)^2) p_arvo21 ## ENNUSTAMINEN ## Otetaan rajatumpi aineisto, joka päättyy vuoden 1981 4. neljännekseen. Toisin sanoen pyritään tekemään "oikea" ennuste. ## Ennusteen lopputuloksessa näkyvät itse ennuste, määritellyn luottamusvälin (ci=0.95) ylä- ja alaraja, sekä CI joka kertoo ## luottamusvälin "leveyden" eli ylä- ja alarajan välisen etäisyyden. Ennustevälit esitetään graafisesti parilla eri komennolla. ## Koodi laskee ennusteen rekursiivisesti kuten monisteen s. 14 on esitetty (eli ensin lasketaan yhden askeleen ennuste, ## jota hyödynnetään toisen askeleen ennustetta laskettaessa jne.) ger1981 = window(ld_ger,end = c(1981,4)) var2_1981 = VAR(ger1981, p = 2, type = c("const"),season = NULL, exogen = NULL, lag.max = NULL) ennuste = predict(var2_1981, n.ahead=4, ci=0.95) ennuste$fcst plot(ennuste) fanchart(ennuste) ## Ks. http://www.jstatsoft.org/v27/i04/paper, s. 12.