Thursday, August 23, 2012

Zaman Serisi Tahmini (Forecasting), R, VAR, rpy2

Zaman serisi modellemesinde VAR modelleri bugunlerde en populer olanlari. VAR, vector autoregression kelimelerinden geliyor, yani tek bir zaman serisi degil, birkac tanesine birden, ayni anda modellemeye ugrasiyoruz.

Autregression, zaman serisinin kendi kendisiyle regresyona sokulmasidir; bilindigi gibi regresyon y = f(x) modellemesinde gurultulu elde edilen y ile, x degerleri arasindaki baglantinin bulunmasina yardim eder (eger f(x) lineer ise iyi sonuclar da bulur). Tek boyutlu zaman serisi modellemesi icin soyle bir numara kullanilir, serinin kopyasi alinir, bir geri kaydirilir, x bu kaydilirilmis seri, y esas seri olur, bu ikili regresyona sokulur. Boylece zaman serisinin kendisini ve regresyon mekanizmasi kullanilarak zaman serisi tahmini yapilabilir.

VAR ise bunu cok boyutlu yapar. Her seriyi hem kendisi, hem de diger tum serilerin p kadar gecmis degeri goz onune alinir. Oldukca guclu bir metottur.

Bu alanda unlu isimlerden Sims'i bilmek gerekir, 1980 yilinda yazdigi ve kendi alanini elestirdigi bir makalede makroekonomide yapisal modeller yerine, ciplak veriye bakmak gerektigini, ve bunu yapmak icin her zaman serilerine tek baslarina degil tum diger serilere de baglantilarini goz onune alarak incelemek gerektigini soyler. VAR matematigi buradan cikmistir. Granger ismi de vardir, VAR modellemesi sonrasi serilerin "birbirine ne kadar etki ettigini" hesaplayan "Granger istatistigi" mesela ona aittir.

Isin matematigine sonra daha detayli girebiliriz, simdilik kodlama acisindan ornekleri verelim. Bu alanda R kodculari cok aktif, o yuzden bir R paketi vars kullanacagiz, ve onu Python uzerinden cagiracagiz.

Diyelim ki bir predict-1.csv icinde bir ulkenin GDP ve tuketim verileri (cons) var, 1959-2009 arasi icin (bu oldukca standart bir veri seti). Once R kurulur

sudo apt-get install r-base-dev r-base python-rpy2

Sonra R'ye girilir

> install.packages("vars")

Simdi su R kodu kullanilabilir

library("vars")

file = "predict-1.csv"
a <- read.csv(file, header = TRUE, sep = ",", na.strings="")

impute.med <- function(x) {
    z <- median(x, na.rm = TRUE)
    x[is.na(x)] <- z
    return(x)
}

a2 <- sapply(a, function(x){
    if(is.numeric(x) & any(is.na(x))){
            impute.med(x)
        } else {
            x
        }
    }
)

out <- VAR(a2, p = 2, type = "const")
out.prd <- predict(out, n.ahead = 30, ci = 0.95)


Bunu Python'dan cagirmak icin rpy2 kullaniriz,

import os, sys
import numpy as np
import rpy2.robjects
from datetime import date, timedelta

f = file("predict.R")
code = ''.join(f.readlines())
result = rpy2.robjects.r(code)
res = [['gdp','cons']]
for i in range(30):
    res.append([str(result[0][0][i]),str(result[0][1][i]) ] )

res = np.array(res)
np.savetxt('predict-2.csv',res,delimiter=",",fmt='%s')


Python isledikten sonra sonuc predict-2.csv icinde olacak. Sonuclar 2009 sonrasi 30 sene sonrasi icin gdp ve tuketim rakamlarini tahmin edecek.

Eger pur Python kullanmak isteseydik, scikits statsmodels adinda bir paketi de kullanabilirdik. Bu durumda hic R kodlamasi olmayacak, kurmak icin

https://github.com/statsmodels/statsmodels

Bu kod

import os
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
 
def pad(data):
    bad_indexes = np.isnan(data)
    good_indexes = np.logical_not(bad_indexes)
    good_data = data[good_indexes]
    interpolated = np.interp(bad_indexes.nonzero()[0], good_indexes.nonzero()[0], good_data)
    data[bad_indexes] = interpolated
    return data

data = np.genfromtxt("predict-1.csv", skip_header=1, delimiter=',')
data = np.apply_along_axis(pad, 0, data)
model = VAR(data)
res = model.fit(2)

f = res.forecast(data[-2:], 30)
np.savetxt('predict-3.csv',f,delimiter=",",fmt='%s')


Ustteki kod sonuclar predict-3.csv icine yazar.

VAR ile zaman serisi tahminlerinde onemli bazi konular incelenen verinin (zaman serisinin) duragan (stationary), ve beraber entegre (co-integrated) olup olmadigidir -- bu durumlarda bazi ek numaralar kullanmak gerekebilir, mesela duragan bir veri seti yoksa serinin farklarini kullanmak gibi..

Ana veri

gdp,cons
2710.349,1707.4
2778.801,1733.7
2775.488,1751.8
2785.204,1753.7
2847.699,1770.5
2834.390,1792.9
2839.022,1785.8
2802.616,1788.2
2819.264,1787.7
2872.005,1814.3
2918.419,1823.1
2977.830,1859.6
3031.241,1879.4
3064.709,1902.5
3093.047,1917.9
3100.563,1945.1
3141.087,1958.2
3180.447,1976.9
3240.332,2003.8
3264.967,2020.6
3338.246,2060.5
3376.587,2096.7
3422.469,2135.2
3431.957,2141.2
3516.251,2188.8
3563.960,2213.0
3636.285,2251.0
3724.014,2314.3
3815.423,2348.5
3828.124,2354.5
3853.301,2381.5
3884.520,2391.4
3918.740,2405.3
3919.556,2438.1
3950.826,2450.6
3980.970,2465.7
4063.013,2524.6
4131.998,2563.3
4160.267,2611.5
4178.293,2623.5
4244.100,2652.9
4256.460,2669.8
4283.378,2682.7
4263.261,2704.1
4256.573,2720.7
4264.289,2733.2
4302.259,2757.1
4256.637,2749.6
4374.016,2802.2
4398.829,2827.9
4433.943,2850.4
4446.264,2897.8
4525.769,2936.5
4633.101,2992.6
4677.503,3038.8
4754.546,3110.1
4876.166,3167.0
4932.571,3165.4
4906.252,3176.7
4953.050,3167.4
4909.617,3139.7
4922.188,3150.6
4873.520,3163.6
4854.340,3117.3
4795.295,3143.4
4831.942,3195.8
4913.328,3241.4
4977.511,3275.7
5090.663,3341.2
5128.947,3371.8
5154.072,3407.5
5191.499,3451.8
5251.762,3491.3
5356.131,3510.6
5451.921,3544.1
5450.793,3597.5
5469.405,3618.5
5684.569,3695.9
5740.300,3711.4
5816.222,3741.3
5825.949,3760.2
5831.418,3758.0
5873.335,3794.9
5889.495,3805.0
5908.467,3798.4
5787.373,3712.2
5776.617,3752.0
5883.460,3802.0
6005.717,3822.8
5957.795,3822.8
6030.184,3838.3
5955.062,3809.3
5857.333,3833.9
5889.074,3847.7
5866.370,3877.2
5871.001,3947.9
5944.020,3986.6
6077.619,4065.7
6197.468,4137.6
6325.574,4203.2
6448.264,4239.2
6559.594,4299.9
6623.343,4333.0
6677.264,4390.1
6740.275,4464.6
6797.344,4505.2
6903.523,4590.8
6955.918,4600.9
7022.757,4639.3
7050.969,4688.7
7118.950,4770.7
7153.359,4799.4
7193.019,4792.1
7269.510,4856.3
7332.558,4910.4
7458.022,4922.2
7496.600,5004.4
7592.881,5040.8
7632.082,5080.6
7733.991,5140.4
7806.603,5159.3
7865.016,5182.4
7927.393,5236.1
7944.697,5261.7
8027.693,5303.3
8059.598,5320.8
8059.476,5341.0
7988.864,5299.5
7950.164,5284.4
8003.822,5324.7
8037.538,5345.0
8069.046,5342.6
8157.616,5434.5
8244.294,5466.7
8329.361,5527.1
8417.016,5594.6
8432.485,5617.2
8486.435,5671.1
8531.108,5732.7
8643.769,5783.7
8727.919,5848.1
8847.303,5891.5
8904.289,5938.7
9003.180,5997.3
9025.267,6004.3
9044.668,6053.5
9120.684,6107.6
9184.275,6150.6
9247.188,6206.9
9407.052,6277.1
9488.879,6314.6
9592.458,6366.1
9666.235,6430.2
9809.551,6456.2
9932.672,6566.0
10008.874,6641.1
10103.425,6707.2
10194.277,6822.6
10328.787,6913.1
10507.575,7019.1
10601.179,7088.3
10684.049,7199.9
10819.914,7286.4
11014.254,7389.2
11043.044,7501.3
11258.454,7571.8
11267.867,7645.9
11334.544,7713.5
11297.171,7744.3
11371.251,7773.5
11340.075,7807.7
11380.128,7930.0
11477.868,7957.3
11538.770,7997.8
11596.430,8052.0
11598.824,8080.6
11645.819,8122.3
11738.706,8197.8
11935.461,8312.1
12042.817,8358.0
12127.623,8437.6
12213.818,8483.2
12303.533,8555.8
12410.282,8654.2
12534.113,8719.0
12587.535,8802.9
12683.153,8865.6
12748.699,8888.5
12915.938,8986.6
12962.462,9035.0
12965.916,9090.7
13060.679,9181.6
13099.901,9265.1
13203.977,9291.5
13321.109,9335.6
13391.249,9363.6
13366.865,9349.6
13415.266,9351.0
13324.600,9267.7
13141.920,9195.3
12925.410,9209.2
12901.504,9189.0
12990.341,9256.0




No comments: