## EPÄSTATIONAARISET AIKASARJAT, kevät 2012
## Koodiesimerkkejä

## 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)


## SIMULOIDAAN SATUNNAISKULKU ESIMERKKINÄ I(1)-PROSESSISTA 

x = rnorm(1000)  ## Innovaation simulointi
z = cumsum(x)  ## Satunnaiskulun muodostaminen (ei vakioita)
plot(z,type = "l")
acf(z) ## Voit muuttaa viipeiden määräksi esimerkiksi 100 kirjoittamalla acf(z, lag=100) 


## YKSIKKÖJUURITESTI

## R:ssä on useita eri vaihtoehtoja
## vars-paketin sisältämä testi

?ur.df 

## ur.df(z, type = c("none", "drift", "trend"), lags = 1, selectlags = c("Fixed", "AIC", "BIC"))
## Sekä viipeiden määrä (”1”) että funktiomuoto ("ei vakiota eikä trendiä", ”vakio” tai (lineaarinen) ”trendi”) valittava.
## Jos käytetään AIC:tä tai BIC:tä, ilmoittaa viipeiden määrä asteen valinnassa käytetyn maksimiasteen.
## Tulostuksen lopussa listatuista testisuureista ”tau-testisuureet” vastaavat kurssilla käsiteltyjä ”t-testisuureita” (niiden arvot näkyvät myös estimointitulosten yhteydessä).
## Phi-testisuureet testaavat ovatko myös vakion tai vakion ja aikatrendin kertoimet nollia (näissäkin tarvittavat jakaumat poikkeavat tavanomaisesta).

## Miten käy satunnaiskulun tapauksessa?

summary(ur.df(z, type = c("none"), lags = 1)) 

## -> nollahypoteesi jäi voimaan


## Pienillä ja pienehköillä havaintomäärillä yksikköjuuritestit eivät ole kovin voimakkaita.
## Havainnollistetaan tätä simuloidulla realisaatiolla AR(1)-mallista.


T = 100 ## Havaintojen lukumäärä
N = 100 ## Toistojen lukumäärä
adf.mat = numeric(N) ## Alustetaan matriisi, johon tallennetaan testisuureen arvoja

Phi_1 = 0.95 ## AR(1)-kerroin
y = epsilon = matrix(,T,N) ## Alustetaan tarvittavat matriisit
for(n in 1:N){
epsilon[,n] = rnorm(T) ## Poimitaan T havaintoa standardinormaalijakaumasta
y[1,n] = epsilon[1,n]

## Simuloidaan T-havainnon AR(1)-malli

for(i in 2:T){
y[i,n] = Phi_1*y[(i-1),n]+epsilon[i,n]
}

adf = ur.df(y[,n], type = c("none"), lags = 1)
adf.mat[n] = adf@teststat[1,1]
}

## 5%:n merkitsevyystason raja on -1.95, kun mallissa ei ole vakiota eikä aikatrendiä
## Montako näistä N:n sarjan yksikköjuuritestistä hylättäisiin 5% tasolla?

length(which(adf.mat < (-1.95)))

## Mitä tulos kertoo testin voimasta? Mitä tapahtuu, jos kasvattaa T:tä? 
## Jos haluat tarkempia tuloksia, voit kasvattaa N:n arvoa (tällöin laskenta kestää pidempään).


#######################################################################


## 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 tiedoston

## 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 1980 I - 2000 IV (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")

##Valitaan tarkasteluun esimerkiksi työttömyyssarja

U=Canada[,"U"]
plot(U,nc=1,xlab="") ## Aikasarjan U kuva 
summary(U) ## Aikasarjasta U laskettuja tunnuslukuja

## Jos haluat kuvat ja tunnusluvut kaikista aineiston aikasarjoista, korvaa edellä ”U” sanalla ”Canada”


## MALLIN ASTEEEN VALINTA

## (1) Mallin asteen valinta auto- ja osittaisautokorrelaatiofunktion avulla. 
## Huomaa, että silloin, kun aikasarjassa on selvä trendi, autokorrelaatiofunktio ei ole informatiivinen. 

corU=acf(U, lag=40, type = "correlation",xlab="Viive", main="Empiirinen autokorrelaatiofunktio työttömyyssarjalle")

pcorU=pacf(U, lag=40, main="Empiirinen osittaisautokorrelaatiofunktio työttömyyssarjalle")

## (2) Mallin asteen valinta mallinvalintakriteerien avulla. 
## select=VARselect(U, lag.max = 6, type = c("const", "trend", "both", "none"), season = NULL, exogen = NULL)
## Näin muotoiltuna valinta kohdistuu yhteen sarjaan eli tässä U:hun. Jos U:n paikalla on Canada, valitaan VAR-mallin aste. Tähän palataan myöhemmin.
## type = c("const", "trend", "both", "none") valitsee mallin tyypin. Valitaan vakiollinen malli:

select=VARselect(U, lag.max = 6, type = c("const"), 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”). 
## Huomaa, että tässä mallin tyypin määrittely poikkeaa hieman aikaisemmasta vastaavasta. 
## Mallissa on lisäksi on mahdollista käyttää kausi-indikaattoreita ja eksogeenisia muuttujia (eksogeenisia muuttujia sisältäviä malleja ei kurssilla käsitellä).


## Huomataan, että sekä autokorrelaatiofunktioiden että mallinvalintakriteerien perusteella päädytään AR(2)-malliin

## YKSKKÖJUUREN TESTAUS

## Testataan nyt yksikköjuurta edellä valitussa mallissa käyttäen vars-paketin testiä. 
## ur.df(y, type = c("none", "drift", "trend"), lags = 1, selectlags = c("Fixed", "AIC", "BIC")) 
## Tulostuksen lopussa listatuista testisuureista ”tau-testisuureet” vastaavat kurssilla käsiteltyjä ”t-testisuureita” (niiden arvot näkyvät myös estimointitulosten yhteydessä).
## Phi-testisuureet testaavat ovatko myös vakion tai vakion ja aikatrendin kertoimet nollia (näissäkin tarvittavat jakaumat poikkeavat tavanomaisesta).

## Viipeiden määrä voidaan siis valita "käsin"

summary(ur.df(U, type = c("drift"), lags = 1))

## Voit kokeilla tässä myös eri asteita ja verrata korkeimman asteen estimaattia keskivirheeseensä tavanomaiseen tapaan. 
## Koska tässä testaus liittyy viivästettyjen differenssien kertoimiin, voidaan menettely osoittaa (asymptoottisesti) päteväksi.
## Huomaa kuitenkin, että yksikköjuurta testatessasi tilanne on toinen, joten (viivästetyn) differensoimattoman muuttujan kohdalla
## tämä menettely ei ole pätevä eikä tulosteen "merkitsevyystähtiin" (*) pidä kiinnittää huomiota (niissä p-arvo lasketaan t-jakaumasta, 
## mikä on tässä tapauksessa väärin; ks. esim HT 3.3 ja 3.4)  

## Viipeiden määrä voidaan valita informaatiokriteerin avulla myös tässä yhteydessä

summary(ur.df(U, type = c("drift"), lags = 6, selectlags=c("BIC")))


## Jos havaitussa aikasarjassa havaitaan yksikköjuuri, voidaan halutessa tutkia onko kysymyksessä mahdollisesti I(2)-prosessin tuottama sarja 
## testaamalla yksikköjuurta differensoidussa  aikasarjassa. Differenssi voidaan muodostaa komennolla

dy = diff(y)