## EPÄSTATIONAARISET AIKASARJAT 2012

## R-koodia sekä stationaariselle että epästationaariselle VAR-MALLILLE; ks.lisää: http://www.jstatsoft.org/v27/i04/paper
## Jossakin kohtaa saatetaan hyödyntää syksyn 2011 moniulotteisten aikasarjojen kurssilla käsiteltyjä koodeja

## TARVITTAVIA PAKETTEJA

install.packages("vars") ## VAR-mallin estimoinnissa
library(vars)
install.packages("tseries") 
library(tseries)
install.packages("fUnitRoots") 
library(fUnitRoots)
library(stats)
library(nlme)


## AINEISTON HAKU (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 tiedosto

## Vaihtoehtoinen menetelmä
Data = read.table("C:/Opiskelu/Tilastotiede/moniulotteiset 2011/Laskuharjoitukset/XXX.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ää. 
Canada = ts(Data[,],start = c(1980,1),deltat = 1/4)

## Käytetään esimerkkinä Kanadasta saatua neljännesvuosittaista talousaineistoa ajalta 1980I-2000IV (sama esimerkki on ylempänä olevan linkin manuaalissa)
## Canadian data: e=log of employment, U=enemplyment rate, rw=log of real wage, prod=log difference between GDP (gross domestic product) and employment
## Aineisto löytyy valmiina R-ohjelmasta, joten edellä mainittua aineiston nimien laittamista tietokantaan ja aikasarjaksi muuntamista ei tässä tarvita 
## (huolehdi, että olet ladannut tarvittavat paketit)

data("Canada")
plot(Canada,nc = 1,xlab = "") ## Aineiston aikasarjojen kuvat 
summary(Canada) ## Aineiston aikasarjoista laskettuja tunnuslukuja

## MALLIN ASTEEN VALINTA MALLINVALINTAKRITEESIEN AVULLA

## Kriteerit ovat Akaiken AIC, Hannanin ja Quinnin HQ, Schwarzin SC=BIC ja Akaiken FPE (Final Prediction Error, jota ei ole käsitelty kurssilla).
## select = VARselect(Canada, lag.max = 6, type = c("const", "trend", "both", "none"), season = NULL, exogen = NULL)
## type = c("const", "trend", "both", "none") valitsee mallin tyypin. Valitaan vakion ja trendin (eli lineaarisen trendin) sisältävä malli:

select = VARselect(Canada, lag.max = 6, type = c("both"), season = NULL, exogen = NULL)
select$selection ## Tulostaa eri kriteerien ehdottamat asteet.
select ## Tulostaa edellisten lisäksi kriteerifunktioiden arvot kaikilla tarkasteltavilla viipeillä.

## Edellä annetaan maksimaalinen viive (”6”), mallin tyyppi (”vakio”, ”trendi”, ”molemmat”, ”ei kumpaakaan”). 
## Mallissa on lisäksi on mahdollista käyttää kausi-indikaattoreita ja eksogeenisia muuttujia (eksogeenisia muuttujia sisältäviä malleja ei kurssilla käsitellä).
## Seuraavassa asteeksi on valittu 2.

## Estimoidaan malli ja tutkitaan residuaaleja myöhempää yi-asteen testausta silmällä pitäen. Tässä malli estimoidaan kuten syksyn 2011 moniulotteisten aikasarjojen 
## kurssilla käyttäen monisteen yhtälön (4.1) esitystä mallista (mahdollisesti aikatrendillä täydennettynä).
## Kun yhteisintegroituvuutta epäillään, ei näin saatuja estimaatteja ole helppo tulkita. Poikkeuksena on korkeimman asteen kerroinmatriisin estimaatti 
## (ks. monisteen s. 50 yhtälöt (5.4) ja (5.5)). Tähän liittyen voidaan syksyn 2011 moniulotteisten aikasarjojen kurssilla esitettyä peräkkäisiin 
## uskottavuusosamäärätesteihin perustuvaa mallin asteen valintamenettelyä myös käyttää sellaisenaan khi^2-jakauman kanssa (ks. monisteen s. 31 loppupuoli ja s. 32 esimerkki).

var_lag2 = VAR(Canada, p = 2, type = c("both"),season = NULL, exogen = NULL, lag.max = NULL) ## Estimointi
plot(var_lag2) ## Residuaaleihin liittyviä kuvia 

## Lisää residuaalikuvia saadaan seuraavasti

res = resid(var_lag2) 
colnames(res) = c("res_e","res_U","res_rw","res_prod")

## Auto- ja ristikorrelaatiot 
acf(res,lag.max = 12, type = "correlation")

## Vastaavat neliöidyille residuaaleille
colnames(res) = c("res_e^2","res_U^2","res_rw^2","res_prod^2")
acf(res^2,lag.max = 12, type = "correlation")  

## Kvantiili-kvantiili kuviot kuvaavat kuinka hyvin residuaalit vastaavat normaalijakauman kertymäfunktion kvantiilipisteitä (yhtälöittäin)
## Mitä enemmän kuvion pisteet poikkeavat 45 asteen suorasta sitä enemmän virheiden komponenttien voidaan epäillä poikkeavan normaalijakaumasta.

qqnorm(res[,1])  ## Valitaan ensimmäisen yhtälön residuaali

## Syksyn 2011 käytetty Portmanteau-testi esi sovellu tässä (sellaisenaan), jos prosessi on I(1).

## YHTEISINTEGROITUVUUSASTEEN TESTAUS 

?ca.jo

## ca.jo(x, type = c("eigen", "trace"), ecdet = c("none", "const", "trend"), K = 2,spec=c("longrun", "transitory"), season = NULL, dumvar = NULL)
## type valitsee käytettävän testin, joka on joko jälkitesti (”trace”) tai maksimiominaisarvotesti (”eigen”).
## ecdet: Kun valitsee "trend", niin mallissa on vakio ja aikatrendi (aikatrendi yi-relaatiossa). Kun valitsee "const", niin mallissa on vakio (yi-relaatiossa). 
## Valittaessa ”none”, käytetään mallia, jossa ei ole vakiota eikä aikatrendiä. Huomaa, että monisteen mallia (5.14) ei ole mahdollista käyttää.
## K on merkintä mallin asteelle eli K on monisteen p.
## spec: Valitse "transitory", niin käytät samaa parametrointia kuin kurssilla (ohjelmassa on mahdollista käyttää myös vaihtoehtoista parametrointia).

## Valitaan asteeksi kaksi ja käytetään jälkitestiä

vecm = (ca.jo(Canada[,c("rw", "prod","e","U")], type = c("trace"), ecdet = c("trend"), K = 2,spec = c("transitory"), season = NULL, dumvar = NULL))
summary(vecm) ## Tämä tulostaa testisuureiden arvot kuten monisteen s. 72 taulukossa (tosin päinvastaisessa järjestyksessä).

## Lisäksi tulostetaan testisuureiden laskemisessa käytettävät (yleistetyt) ominaisarvot ja ominaisvektorit.
## Jälkimmäiset ohjelma normalisoi automaattisesti ensimmäisen valitun sarakkeen suhteen (tässä siis "rw":n suhteen).
## Testisuure hylkää oletuksen nollasta yi-asteesta mutta jättää voimaan asteen yksi.

## VIRHEENKORJAUSMALLIN ESTIMOINTI VALITULLA YI-ASTEELLA

## Kun edellisten perusteella on valittu yi-aste, niin seuraavilla komennoilla voidaan estimoida valittu virheenkorjausmalli
## (yhteisintegroitunut VAR) pns-menetelmällä (eli monisteen s. 76 puolenvälin jälkeen esitettyllä pns-menetelmällä).

?cajorls

## r on yi-aste
vecm.r1 = cajorls(vecm,r = 1)
summary(vecm.r1$rlm) ## Estimointitulokset yhtälöittäin

## Muodostetaan alpha- ja betametriisien estimaatit ja lasketaan niiden alkoiden keskivirheet ja "t-suhteet" .
## Huomaa, että alphan alkioiden estimaattien keskivirheet ja "t-suhteet" saadaan myös edellisistä tuloksista, 
## mutta ne poikkeavat hieman seuraavista (luultavasti eri laskutavasta johtuen).

alpha = coef(vecm.r1$rlm)[1, ] ## Latausmatriisi alpha
beta = vecm.r1$beta ## Yhteisintegroituvuusmatriisi beta
resid = resid(vecm.r1$rlm) ## Residuaalit
N = nrow(resid)  ## Residuaalisarjojen pituus
omega = crossprod(resid)/N  ## Virhekovarianssimatriisin Omega estimaatti
alpha.se = sqrt(solve(crossprod(cbind(vecm@ZK %*% beta, vecm@Z1)))
                               [1,1]*diag(omega)) 
alpha.t = alpha/alpha.se 
beta.se = sqrt(diag(kronecker(solve(crossprod(vecm@RK[,-1])),
                solve(t(alpha) %*% solve(omega) %*% alpha))))
beta.t = c(NA, beta[-1]/beta.se) 

alpha ## alphan estimaatti
alpha.se ## alphan estimaatin alkioiden keskivirheet
alpha.t ## alphan estimaatin alkioiden "t-suhteet"

beta ## betan estimaatti
beta.se ## betan estimaatin alkioiden keskivirheet
beta.t ## betan estimaatin alkioiden "t-suhteet"

## Residuaalitarkasteluja

resid = ts(resid)
plot(resid) ## Residuaalit graafisesti
acf(resid, lag.max = 12) ## Residuaalien auto- ja ristiorrelaatiot
acf(resid^2, lag.max = 12) ## Neliöityjen residuaalien auto- ja ristiorrelaatiot
qqnorm(resid[,4])  ## Kvantiili-kvantiilikuvio (valitaan ensimmäisen yhtälön residuaali).

## IMPULSSIVASTEET

## Virheenkorjausmalli muunnetaan ensin ”tavalliseen VAR-muotoon”. Tapauksessa, jossa ei ole vakiota tai aikatrendiä tämä vastaa  
## monisteen yhtälön (5.7) muuttamista yhtälön (5.4) muotoon, jolloin yhtälön (5.4) parametrimatriisit siis toteuttavat yi-rajoitteet.

vecm.level = vec2var(vecm,r = 1) ## Virheenkorjausmallin muuntaminen ”tavalliseen VAR-muotoon”

## Jos halutaan ortogonaaliset impussivasteet, niin ortho = TRUE, muuten ortho = FALSE. 
## Argumentti ”n.ahead” ilmoittaa laskettavien impulssivasteiden lukumäärän.
## Luottamusvälit muodostetaan (eräällä) bootstrap-menetelmällä, jonka tuloksia voidaan tarkentaa valitsemalla ”runs” suuremmaksi.

impuls = irf(vecm.level, n.ahead = 15, ortho = FALSE, boot = TRUE,ci = 0.95,runs = 100, cumulative = FALSE)  ## Lasketaan impulssivasteet
plot(impuls)  ## Piirretään impulssivasteiden kuvat