Zpět na hlavní stránku

V rámci své diplomové práce jsem vytvořil aplikaci běžící pod Matlabem, která umožňuje pracovat s různými typy identifikace při práci s modely vektorové autoregrese (VAR). Nejde o plnohodnotnou (uživatelsky přívětivou) aplikaci - spíše se jedná o soubor různých skriptů, které na sebe navazují. Aplikace taktéž neobsahuje žádné prvky GUI.

Aplikaci jsem programoval ve studentské verzi Matlabu z roku 2007 (tuším verze 7.4). Jelikož se práce s VAR modely opírá vesměs o maticovou algebru, předpokládám, že je aplikace spustitelná i v jiných verzích Matlabu.

Pro spuštění aplikace je třeba mít k dispozici Matlab včetně modulu "Symbolic Math Toolbox". Pokud tento modul nemáte nainstalovaný, pak bude plně funkční jen identifikace pomocí Choleskiho dekompozice a následná analýza funkcí odezvy.

OBSAH:

Instalace
Data
Inicializace
Odhad parametrů
Funkce odezvy
Jiný typ identifikace
Identifikace Blanchard-Quah

INSTALACE

Po stažení souboru 'VARtools.zip' je nutno tento soubor dekomprimovat na libovolné místo na pevném disku (např. C:\Matlab\VAR).

Při každém spuštění se v Matlabu musí přenastavit 'Current Directory' tak, aby se cesta shodovala s adresářem, do kterého se toolbox předtím dekomprimoval (čili taktéž C:\Matlab\VAR).

adresar

V tomto adresáři by měl být mimo jiné soubor 'activate.m'. Aby aplikace fungovala správně, je nutné nastavit cestu k adresáři s instalací i v tomto souboru pomocí příkazu 'edit activate' - cesta se přenastaví změnou řetězce na druhém řádku. Ostatní řádky by měly zůstat nezměněny. Takto modifikovaný soubor je třeba uložit a v budoucnu jej není třeba znovu přenastavovat.

V tuto chvíli je aplikace zcela připravena k výpočtům.

nahoru

DATA

Ve své diplomové práci jsem pracoval se čtvrtletními daty za období 1Q1995 až 2Q2007 (úrovňová data, stacionarizovaná data ve formátu .xls). Aplikace obsahuje časové řady pro HDP a M2 v reálném vyjádření, dále kurs CZK/EUR, PRIBOR a cenový index (CPI). Samozřejmě je možné použít jakákoliv vlastní data, ale je třeba pamatovat na to, že aplikace umí pracovat maximálně s pětirovnicovým VAR modelem. Postup zadání vlastních dat je možno nalézt zde.

Příkazem 'edit data' se volí, jaké proměnné budou zařazeny do modelu. Symbol '%' odstraníme u té proměnné, kterou chceme v modelu ponechat a naopak se tento symbol připíše u té proměnné, kterou chceme vynechat.

Např. následující schéma ponechá v modelu pouze HDP a M2, model bude tedy dvourovnicový.

VYRAZENI PROMENNYCH Z MODELU
%1-Y, 2-M2, 3-kurs, 4-ceny, 5-PRIBOR

vec_prom = 1:5;
data_y(:,5) = [];data_levels(:,5) = [];vec_prom(5) = [];
data_y(:,4) = [];data_levels(:,4) = [];vec_prom(4) = [];
data_y(:,3) = [];data_levels(:,3) = [];vec_prom(3) = [];
%data_y(:,2) = [];data_levels(:,2) = [];vec_prom(2) = [];
%data_y(:,1) = [];data_levels(:,1) = [];vec_prom(1) = [];

Poslední částí v přípravě dat je nastavení permutační matice (PERM), která stanoví pořadí proměnných v modelu, což je důležité při identifikaci modelu pomocí Choleskiho dekompozice. Pokud tuto matici nastavíme na jednotkovou, pak v návaznosti na uvažovaný příklad proměnná HDP bude zařazena do modelu jako první, M2 jako druhá, ostatní proměnné jsou ignorovány.

untitled2
%PORADI PROMENNYCH (i-tou pozici zabere j-ta promenna)
%permutacni matice - pouze blok o rozmeru
        %poctu promennych v modelu
        %je bran v uvahu (levy horni blok
        %matice PERM), zbytek permutacni
        %matice je ignorovan,
        %kazdy radek a sloupec obsahuje prave jednu '1'
PERM = [1 0 0 0 0;
        0 1 0 0 0;
        0 0 1 0 0;
        0 0 0 1 0;
        0 0 0 0 1];

Pokud chceme pořadí proměnných změnit, stačí v permutační matici editovat jen levý horní blok 2 krát 2, který přísluší uvažovanému dvourovnicovému modelu s proměnnými HDP a M2. Ostatní řádky a sloupce matice PERM jsou ignorovány:

untitled2
PERM = [0 1 0 0 0;
        1 0 0 0 0;
        0 0 1 0 0;
        0 0 0 1 0;
        0 0 0 0 1];

nahoru

INICIALIZACE

Máme-li v modelu příslušná data, je třeba nastavit celou řadu parametrů modelu, což provedeme nejlépe příkazem 'edit ini'.

ini
konstanta = 1;              %1 = ano, 0 = ne
trend = 1;                  %1 = ano, 0 = ne
sezonni_dummy = 1;          %1 = ano, 0 = ne
strukturalni_zlom = 8;      %0 = zadny zlom, jinak cislo pozorovani, ve kterem dochazi ke zlomu
pouzit_exogenni_prom = 1;   %1 = ano, 0 = ne
lag = 3;                    %max. delka zpozdeni
alfa = 0.1;                 %hladina vyznamnosti
reakci = 20;                %na kolik obdobi impuls-reakce
simulaci = 1500;            %kolikrat simulovat reakce na impuls
odsazeni = 0.04;            %okraj grafu 1 = 100% rozpeti

V tomto skriptu tedy nastavíme, jestli model obsahuje konstantu, trendovou proměnnou, umělé sezónní proměnné, atd.

Parametr 'strukturalni_zlom' lze nastavit na nulu v případě, že předpokládáme stabilitu parametrů modelu v celém sledovaném období. Pokud předpokládáme variabilitu regresních parametrů v čase, je možno tento parametr nastavit na číslo pozorování, ve kterém uvažujeme strukturální zlom. Na datech za ČR jde zajisté o období v letech 1997-1998.

Parametr 'pouzit_exogenni_prom' je možno nastavit na 1, pokud máme k dispozici hodnoty nějaké exogenní proměnné, kterou přirozeně mohou být i pomocné 0-1 proměnné. Data pro tyto proměnné jsou v matici 'data_dummy.mat'.

'Lag' je důležitým parametrem, pomocí kterého se nastaví maximální délka zpoždění. Pokud potřebujeme přihlédnout při volbě tohoto parametru k informačním kritériím, je možné použít příkazu 'aic', který se napíše do dialogového okna Matlabu.

Poslední tři parametry jsou důležité při analýze funkcí odezvy. Jde o délku reakčního období, počet simulačních iterací, pomocí kterých se konstruují konfidenční intervaly a parametr 'odsazení' určí, jak jsou data v grafickém vyjádření vzdálena od okraje grafu (v procentech rozměrů grafu).

Po nastavení všech parametrů modelu je třeba soubor 'ini.m' uložit a následně spustit v dialogovém okně Matlabu příkazem 'ini'.


nahoru

Odhad parametrů

Po provedení inicializace prostřednictvím spuštění skriptu 'ini', je možné odhadnout regresní parametry a doprovodné statistiky pomocí příkazu 'estimate'. Např. odhadnuté parametry lze zobrazit pomocí příkazu 'B_all', p-hodnoty pomocí 'p_value', atp.

Pokud je model špatně specifikován - mimo jiné vzhledem k počtu stupňů volnosti, zobrazí se hláška, že je model nestabilní (např. byl nalezen jednotkový kořen).

Z testů diagnostické kontroly je v aplikaci implementován pouze test normality reziduí, neboť ten v PcGivu je závislý na pořadí proměnných v modelu. Tento test vyvoláme příkazem 'jb'. Grafickou analýzu reziduí pomocí funkcí ACF/PACF je možno vyvolat příkazem 'acf'.

Pomocí funkce 'chow' lze testovat, v jakém pozorování dochází nejspíše ke strukturálnímu zlomu. Tato funkce ale není aplikovatelná na data v tomto modelu, protože obsahují příliš málý počet pozorování.


nahoru

Analýza funkcí odezvy

Standardně je použita při identifikaci modelu SVAR Choleskiho dekompozice, lze však využít i jiný typ identifikace (viz. zde).

Příkaz 'impulse_response' vypočítá na základě regresních parametrů funkce odezvy a zobrazí je ve známé grafické podobě. Výsledek může vypadat v modelu s proměnnými HDP a M2 takto:

funkce odezvy

Chceme-li zobrazit funkce odezvy včetně příslušných konfidenčních intervalů, je třeba model simulovat, což zabere pár sekund až desítek sekund. Simulaci je možné provést pomocí 'simulace_stabilita'. Hladinu významnosti a ostatní parametry potřebné pro simulaci je možno nastavit pomocí 'edit ini'. Ve výstupu simulace je uveden minimální potřebný počet simulačních iterací. Pokud je skutečný počet simulací mnohem menší než toto číslo, je potřeba v inicializaci modelu počet simulací zvýšit a provést celou simulaci znovu, jinak je ve výsledku zkreslena hladina významnosti, na které se činí závěry.

Výsledek může vypadat např. takto:

funkce odezvy

Chceme-li prezentovat pouze hranice konfidenčních intervalů (nikoliv celé plochy), lze použít příkaz 'bands' (ne rovnou, ale až po simulaci!). Funkce odezvy pak vypadají takto:

funkce odezvy

nahoru

Jiný typ identifikace

Při identifikaci obecného SVAR modelu, který není rekurzivního typu, je nutné řešit systém nelineárních rovnic. Ve své aplikaci jsem implementoval algoritmus založený na Newtonově-Raphsonově metodě gradientního spádu, který funguje pro přesně identifikované modely.

Pomocí příkazu 'edit identifikace' je třeba omezit některé koeficienty strukturní matice. V rekurzivním modelu jsou omezeny koeficienty nad (pod) diagonálou matice 'I', obecně mohou být omezení jakákoliv (model musí ovšem zůstat přesně identifikován). Následně je třeba tento soubor uložit!! Identifikace v případě třírovnicového modelu může vypadat např. takto:

ident
%0 = parametr matice B je omezen (nulovy)
%1 = parametr dopocitan v ramci modelu
%   je-li model napr. rozmeru 2,
%   prvky od 3. dimenze nebudou brany v uvahu,
%   neni treba zde prenastavovat tyto prvky na 0
%je treba, aby model byl presne identifikovan!!!

I = [1 0 1 0 0;
     0 1 1 0 0;
     1 0 1 1 0;
     0 0 1 1 1;
     0 1 1 1 1];

max_iteraci = 40;
%maximum iteraci vypoctu (staci male cislo)
%lepe nastavit mene a pokud presnost neni
%dostatecna teprve zvysovat (presnost znamena, ze sigma_strukt
%je na vystupu +/- diagonalni
%pokud je vypocet prilis dlouhy, znaci to neefektivnost algoritmu
%bud zkusit znovu s jinymi random cisly (vetsinou nepomuze),
%nebo zmenit identifikaci (bohuzel)

Podobně jako u výše zmíněné permutační matice, i v matici I se ignorují ty řádky a sloupce, které jsou v modelu navíc. Je-li model třírovnicový, potom se 4. a 5. řádky a sloupce matice I nemusí vůbec přenastavovat (stejně je nakonec skipt nebere v úvahu).

Čím více se strukturní omezení v matici I neslučují s rekurzivním modelem (pokud se omezení matice I výrazně liší od trojúhelníkové struktury), tím déle trvá výpočet strukturních parametrů. Navíc algoritmus, který jsem implementoval do aplikace, nepatří mezi nejlepší, proto se často stane, že za výsledek je označeno řešení, kdy kovarianční matice strukturního tvaru není diagonální. V takovém případě je možné zvýšit počet numerických iterací (parametr 'max_iteraci' v panelu identifikace pod maticí I), nebo model identifikovat nějak jinak (bohužel).

To, zda je identifikace přijatelná poznáme tak, že spustíme 'svar' a ve výstupu by měla matice 'sigma_strukt' být přibližně diagonální.

Následně je možná analýza funkcí odezvy pomocí výše zmíněných příkazů - 'impulse_response', 'simulace_stabilita' a 'bands'.


nahoru

Identifikace Blanchard-Quah

Skript pro identifikaci typu BQ není plně automatizován. Sekvence příkazů je násleující: 'ini', 'estimate', 'bq_ini'. V okně Matlabu by v této fázi měly být spočítané matice 'mat_soucet' a 'sigma'. Následně pomocí 'edit fun_bq' je třeba do modelu zadat dlouhodobé omezení.

V tomto skriptu je předvoleno omezení pro dvou a třírovnicový model (u jednoho z modelů je nutné připsat na začátky řádků symbol '%', aby je Matlab ignoroval). Na základě spočtených hodnot matic 'mat_soucet' a 'sigma' ve výstupu okna Matlabu je potřeba do skriptu přímo vepsat rozptyly a kovariance (z matice 'sigma') a dále některé prvky matice 'mat_soucet' - komentář přímo uvnitř skriptu 'fun_bq' +/- vysvětluje které prvky.

V následujícím obrázku je správně vyplněný skript 'fun_bq' (vlevo) na základě hodnot spočtených v matlabovském okně (vpravo), v tomto případě bylo potřeba doplnit 2 hodnoty rozptylů a 1 kovarianci z matice 'sigma' a dvojici prvků matice 'mat_soucet':

Blanchard-Quah
nahoru