Bu blog yazısı Eview blogdan Türkçeye çevrilmiştir ve kodlar R’a çevrilmiştir

Posted by IHSEViews

Giriş

Otoregresif Dağıtılmış Gecikme (ARDL) modelleri, hem bağımlı hem de bağımsız değişkenlerin sadece eşzamanlı olarak değil, aynı zamanda geçmiş (gecikmeli) değerleriyle de ilişkili olduğu doğrusal zaman serisi modelleridir. Özellikle, \(y_t\) bağımlı değişken ve \(x_1, \ldots, x_k\) \(k\) adet açıklayıcı değişken olmak üzere, genel bir ARDL\((p,q_1,\ldots,q_k)\) modeli şu şekilde ifade edilir:

\[ y_t = a_0 + a_1t + \sum_{i=1}^p{\psi_i y_{t-i}} + \sum_{j=1}^{k}\sum_{l_j=0}^{q_j}{\beta_{j,l_j}x_{j,t-l_j}} + \epsilon_t \]

Burada, \(\epsilon_t\) genellikle inovasyonları temsil eder, \(a_0\) sabit terimdir ve \(a_1\), \(\psi_i\) ve \(\beta_{j,l_j}\) sırasıyla lineer trend, \(y_t\)’nin gecikmeleri ve \(k\) adet \(x_{j,t}\) regresörünün gecikmeleri ile ilişkili katsayılardır.

Alternatif olarak, \(L\) gecikme operatörünü belirtmek için kullanılırsa ve \(\psi(L)\) ve \(\beta_j(L)\) gecikme polinomları olarak tanımlanırsa:

\[ \psi(L) = 1 - \sum_{i=1}^{p}{\psi_iL^i} \quad \text{ve} \quad \beta_j(L) = \sum_{l_j=0}^{q_j}{\beta_{j,l_j}L^{l_j}} \]

O zaman, yukarıdaki denklem şu şekilde de yazılabilir:

\[ \psi(L)y_t = a_0 + a_1t + \sum_{j=1}^{k}{\beta_{j}(L)x_{j,t}} + \epsilon_t \]

ARDL modelleri, ekonometrik analizlerde onlarca yıldır kullanılmaktadır, ancak son yıllarda eşbütünleşme ilişkilerini incelemede popülerlik kazanmıştır. Bu bağlamda, Pesaran ve Shin (1998) ile Pesaran, Shin ve Smith (2001) önemli katkılar yapmıştır. Özellikle, ARDL modellerinin, ilgili değişkenlerin entegrasyon derecelerinin yanlış belirtilmesine karşı doğuştan gelen bir sağlamlığa sahip olduğu ve bu nedenle eşbütünleşme analizinde avantajlı olduğu belirtilmiştir.

Bu bağlamda üç durum dikkate alınabilir:

Tüm değişkenler \(I(d)\) (durağanlık derecesi) ve eşbütünleşik değildir. Bu durumda, \(d=0\) ise seviyelerde, \(d>0\) ise uygun farklarda en küçük kareler tekniği kullanılabilir.
Tüm değişkenler \(I(1)\) ve eşbütünleşiktir. Bu durumda:
- \(y_t\)’yi \(x_{j,t}\)’lere karşı seviyelerde regresyon yaparak eşbütünleşik (uzun dönem) ilişki tahmin edilebilir.
- Uygun hata düzeltme modeli (ECM) kullanılarak kısa dönem dinamiklerinin eşbütünleşik ilişkiye uyum hızı tahmin edilebilir.
Bazı değişkenler \(I(0)\), diğerleri \(I(1)\) ve bunlar arasında bazıları eşbütünleşiktir.

Özellikle üçüncü durumda, geleneksel eşbütünleşme yöntemleri genellikle başarısız olur çünkü tüm değişkenlerin aynı entegrasyon derecesine sahip olması gerekir. Bu, her bir değişken için birim kök testlerinin yapılmasını gerektirir ki bu testler, boyut ve güç sorunlarına duyarlıdır.

Alternatif olarak, PSS(2001) sınır testi, bu sınırlamalara tabi değildir ve üçüncü durumun nüanslarını kolayca kapsar. Test, temel olarak, altta yatan vektör otoregresyon (VAR) modelinin ECM’sindeki uzun dönemli değişkenler üzerindeki bir parametre anlamlılık testidir ve tüm veya bazı değişkenler \(I(0)\), \(I(1)\) veya hatta karşılıklı olarak eşbütünleşik olduğunda çalışır. Bir ECM ile bir ARDL modeli arasında bire bir bir ilişki olduğundan ve ARDL modelleri tanıdık en küçük kareler teknikleri kullanılarak tahmin edilip yorumlandığından, ARDL modelleri de facto tercih edilen yöntem haline gelmiştir.

Modelin Belirlenmesi (Specification)

Bir ARDL modeli belirlenirken, bağımlı değişkenin gecikmeleri ve açıklayıcı değişkenlerin gecikmeleri dikkate alınarak bir yapı oluşturulur. Bu model, değişkenler arasındaki kısa ve uzun dönem ilişkileri incelemek için kullanılır. Genel bir ARDL\((p,q_1,\ldots,q_k)\) modeli şu şekilde ifade edilir:

\[ y_t = a_0 + a_1t + \sum_{i=1}^{p}{\psi_i y_{t-i}} + \sum_{j=1}^{k}\sum_{l_j=0}^{q_j}{\beta_{j,l_j}x_{j,t-l_j}} + \epsilon_t \]

Burada: - \(y_t\) bağımlı değişken, - \(x_{j,t}\) bağımsız değişkenler, - \(p\) bağımlı değişkenin gecikme uzunluğu, - \(q_j\) her bir bağımsız değişken için belirlenen gecikme uzunluğudur.

Bu model, değişkenlerin zaman içindeki etkileşimlerini ve uzun dönem ilişkilerini analiz etmeye olanak tanır.

Genel ARDL modeli bu denklem ile tanımlanmış olsa da, üç alternatif temsili mevcuttur. Bu üç temsilin tümü parametre tahmini için kullanılabilir. Ancak:

Birinci temsil, genellikle zamana yayılmış dinamik tahmin için kullanılır.
İkinci temsil, tahmin sonrası uzun dönem (denge) ilişkisinin türetilmesi için kullanılır.
Üçüncü temsil, denklem’in PSS (2001) sınır testi bağlamında koşullu hata düzeltme (CEC) gösterimine indirgenmiş halidir; bkz. Banerjee ve diğerleri (1993).

Bu üç temsilin tümü, bazı öncül sonuçları gerektirir. Ünlü Beveridge-Nelson ayrıştırma prensiplerine dayanarak, \(\psi(L)\) ve \(\beta_j(L)\)’nin her zaman aşağıdaki gibi ayrıştırılabileceğini hatırlayalım:

\[ \psi(L) = \psi(1) + (1 - L) \tilde{\psi}(L) \quad \text{and} \quad \beta_j(L) = \beta_j(1) + (1 - L) \tilde{\beta}_j(L) \]

Burada,

\[ \tilde{\psi}(L) = \sum_{i=0}^{p-1} \tilde{\psi}_i L^i, \quad \tilde{\psi}_i = - \sum_{r=i+1}^{p} \psi_r \]

\[ \tilde{\beta}_j(L) = \sum_{l_j=0}^{q_j-1} \tilde{\beta}_{j,l_j} L^{l_j}, \quad \tilde{\beta}_{j,l_j} = - \sum_{s=l_j+1}^{q_j} \beta_{j,s} \]

\[ \psi(1) = 1 - \sum_{i=1}^{p} \psi_i, \quad \beta_j(1) = \sum_{l_j=0}^{q_j} \beta_{j,l_j} \]

Ayrıca, aşağıdaki bağıntıyı göz önünde bulunduralım:

\[ \psi(L) = 1 - \psi^*(L) \quad \text{burada} \quad \psi^*(L) = \sum_{i=1}^{p} \psi_i L^i. \]

Bunun yanı sıra:

\[ \psi^*(L) = \sum_{i=1}^{p} \psi_i L^i = \left( \sum_{i=1}^{p} \psi_i L^{i-1} \right) L = \left( \psi^*(1) + (1 - L) \tilde{\psi}^*(L) \right) L \]

Burada,

\[ \tilde{\psi}^*(L) = \sum_{i=1}^{p-1} \tilde{\psi}^*_i L^{i-1}, \quad \tilde{\psi}^*_i = - \sum_{r=i+1}^{p} \psi_r \]

\[ \psi^*(1) = \sum_{i=1}^{p} \psi_i. \]

Son olarak, herhangi bir \(z_t\) serisi için her zaman aşağıdaki gösterim yazılabilir:

\[ z_t = z_{t-1} + \Delta z_t \]

İlk Temsil: Zamana Yayılmış Dinamik Regresyon (Intertemporal Dynamics Regression)

ARDL modelinin ilk temsili, bağımlı değişkenin kendi gecikmelerini ve bağımsız değişkenlerin gecikmelerini içeren standart bir regresyon modelidir. Aşağıdaki gibi yazılabilir:

\[ y_t = a_0 + \sum_{i=1}^{p}{\psi_i y_{t-i}} + \sum_{j=1}^{k} \sum_{l_j=0}^{q_j} \beta_{j,l_j} x_{j,t-l_j} + \epsilon_t \]

Bu form, zaman içindeki kısa dönem dinamiklerini açıkça gösterir. Ancak, uzun dönem ilişkileri ve hata düzeltme mekanizmalarını doğrudan içermez.

ARDL uygulamalarının çoğu için tipik başlangıç noktası, zamansal dinamiklerin (intertemporal dynamics) tahmin edilmesidir. Bu formda, \(y_t\) değişkeninin kendi gecikmeli değerleriyle ve aynı zamanda \(k\) adet açıklayıcı değişkenin \(x_{j,t}\) eşzamanlı ve gecikmeli değerleriyle olan ilişkisini tahmin etmekle ilgilenilir. Aslında bu, PS (1998) tarafından incelenen ARDL modelinin temelini oluşturur. Özellikle, denklemi şu şekilde ifade ederiz:

\[ y_t = a_0 + a_1 t + \sum_{i=1}^{p} \psi_i y_{t-i} + \sum_{j=1}^{k} \beta_j(L)x_{j,t} + \epsilon_t \]

\[ = a_0 + a_1 t + \sum_{i=1}^{p} \psi_i y_{t-i} + \sum_{j=1}^{k} \left( \beta_j(1) + (1 - L) \tilde{\beta}_j(L) \right) x_{j,t} + \epsilon_t \]

\[ = a_0 + a_1 t + \sum_{i=1}^{p} \psi_i y_{t-i} + \sum_{j=1}^{k} \beta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\beta}_j(L) \Delta x_{j,t} + \epsilon_t \]

Burada birinci fark operatörü \(\Delta = (1 - L)\) olarak tanımlanmıştır. Denklem doğrudan \(y_t\) için çözülmediğinden, genellikle zamansal dinamiklerin tahmini için bir regresyon modeli olarak yorumlanır. Elbette, yukarıdaki model teorik katsayılar kullanırken, pratik bir regresyon analizi bağlamında aşağıdaki gibi ifade edilir:

\[ y_t = a_0 + a_1 t + \sum_{i=1}^{p} b_{0,i} y_{t-i} + \sum_{j=1}^{k} b_j x_{j,t} + \sum_{j=1}^{k} \sum_{l_j=1}^{q_j-1} c_{j,l_j} \Delta x_{j,t-l_j} + \epsilon_t \]

İkinci Temsil: Uzun Dönem Dinamiklerinin Regresyon Sonrası Türetilmesi

ARDL modelinin ikinci temsili, uzun dönem dinamiklerinin türetilmesine odaklanır. Uzun dönem dengesini ifade eden denklemi elde etmek için, modeldeki tüm değişkenlerin zaman içindeki durağan duruma (steady state) ulaşması gerektiği varsayılır. Bu durumda:

\[ y_t = y_{t-1} = y^*, \quad x_{j,t} = x_{j,t-1} = x_j^* \]

olarak kabul edilir ve hata terimi \(\epsilon_t\) sıfır olur. Modelin uzun dönem denge formu şu şekilde yazılabilir:

\[ y^* = \frac{a_0 + \sum_{j=1}^{k} \sum_{l_j=0}^{q_j} \beta_{j,l_j} x_j^*}{1 - \sum_{i=1}^{p} \psi_i} \]

Bu ifade, bağımlı değişkenin uzun dönem dengesini belirleyen temel faktörleri gösterir.

İkinci gösterim, esasen \(y_t\) ile \(k\) adet açıklayıcı değişken arasındaki uzun dönem ilişkisini türetme girişimidir. Bu nedenle, gösterim \(y_t\)’yi \(x_{j,t}\) cinsinden çözer. Özellikle, denklem (1)’den başlayarak şu ifadeyi dikkate alalım:

\[ \psi(L) \Delta y_t = (1 - L) \psi(L) y_t \]

\[ = (1 - L) \left( y_t - \sum_{i=1}^{p} \psi_i y_{t-i} \right) \]

\[ = (1 - L) \left( a_0 + a_1 t + \sum_{j=1}^{k} \beta_j(L) x_{j,t} + \epsilon_t \right) \]

\[ = a_1 + \sum_{j=1}^{k} \beta_j(L) \Delta x_{j,t} + \Delta \epsilon_t \]

Bundan sonra, \(\psi(L)\)’nin terslenebilir olduğunu, yani karakteristik polinomun

\[ 1 - \sum_{i=1}^{p} \psi_i z^i = 0 \]

köklerinin olduğunu ve böylece varsayarsak, şu ifade geçerli olur:

\[ \Delta y_t = \psi^{-1}(L) \left( a_1 + \sum_{j=1}^{k} \beta_j(L) \Delta x_{j,t} + \Delta \epsilon_t \right) \]

Ayrıca, \(\psi(L) y_t = \psi(1) y_t + \tilde{\psi}(L) \Delta y_t\) olduğunu göz önünde bulundurarak, denklem (1)’i şu şekilde yeniden yazabiliriz:

\[ \psi(1)y_t = \psi(L)y_t - \tilde{\psi}(L)\Delta y_t \]

\[\psi(1)y_t = a_0 + a_1 t + \sum_{j=1}^{k} \beta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\beta}_j(L)\Delta x_{j,t} - \tilde{\psi}(L)\Delta y_t + \epsilon_t\]

\[ = a_0 + a_1 t + \sum_{j=1}^{k} \beta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\beta}_j(L)\Delta x_{j,t} - \tilde{\psi}(L)\psi^{-1}(L) \left(a_1 + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \Delta\epsilon_t \right) + \epsilon_t \]

\[ = a_0^* + a_1 t + \sum_{j=1}^{k} \beta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\beta}_j^*(L)\Delta x_{j,t} + \epsilon_t^* \]

where

\[ a_0^* = a_0 - \tilde{\psi}(L)\psi^{-1}(L)a_1 \] \[= a_0 - \tilde{\psi}(1)\psi^{-1}(1)a_1\]

\[\tilde{\beta}_j^*(L) = \tilde{\beta}_j(L) - \tilde{\psi}(L)\psi(L)^{-1}\beta_j(L)\]

\[\epsilon_t^* = \epsilon_t - \tilde{\psi}(L)\psi(L)^{-1}\Delta\epsilon_t\]

At last, the second representation is formulated as:

\[y_t = \psi^{-1}(1) \left( a_0^* + a_1 t + \sum_{j=1}^{k} \beta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\beta}_j^*(L)\Delta x_{j,t} + \epsilon_t^* \right)\]

\[= a_0 + a_1t + \sum_{j=1}^{k} \theta_j(1)x_{j,t} + \sum_{j=1}^{k} \tilde{\theta}_j(L)\Delta x_{j,t} + \xi_t \tag{5}\]

Burada

\[\alpha_0 = \psi^{-1}(1) \left( a_0 - \tilde{\psi}(1)a_1 \right)\] \[\alpha_1 = \psi^{-1}(1)a_1\] \[\theta_j(1) = \psi^{-1}(1)\beta_j(1)\] \[\tilde{\theta}_j(L) = \psi^{-1}(1) \left( \tilde{\beta}_j(L) - \tilde{\psi}(L)\psi(L)^{-1}\beta_j(L) \right)\] \[\xi_t = \psi^{-1}(1) \left( \epsilon_t - \tilde{\psi}(L)\psi(L)^{-1}\Delta\epsilon_t \right)\]

Denklem (3)’ten, genellikle \(j = 1, \ldots, k\) için \(\alpha_1\) ve \(\theta_j(1)\) ile ifade edilen uzun dönem (trend) parametreleriyle ilgileniyoruz. Aslında, (3) ve denklem (5)’te elde edilen parametre tahminleri arasındaki birebir ilişki göz önüne alındığında, regresyon modeli (4) tahmin edildikten sonra, yukarıdaki parametre formüllerini kullanarak uzun dönem parametrelerinin tahminlerini elde etmek mümkündür. Özellikle, eğer \(\hat{a}_1, \hat{b}_{0,1}, \ldots, \hat{b}_{0,p}, \hat{b}_1 \ldots, \hat{b}_k\) regresyon modelinden (4) elde edilen tahmini katsayıların ilgili alt kümesini gösteriyorsa, uzun dönem parametrelerinin regresyon sonrası tahmini aşağıdaki gibi elde edilir:

\[\hat{\alpha}_1 = \frac{\hat{a}_1}{1 - \sum_{i=1}^{p} \hat{b}_{0,i}} \quad \text{ve} \quad \hat{\theta}_j(1) = \frac{\hat{b}_j}{1 - \sum_{i=1}^{p} \hat{b}_{0,i}}\]

Üçüncü Temsil: Koşullu Hata Düzeltme Formu ve Sınır Testi

ARDL modelinin üçüncü temsili, hata düzeltme mekanizmasını içerir ve değişkenler arasında uzun dönem eşbütünleşme ilişkisi olup olmadığını belirlemek için kullanılır. Bu form şu şekilde yazılabilir:

\[ \Delta y_t = a_0 + \sum_{i=1}^{p-1} \psi_i \Delta y_{t-i} + \sum_{j=1}^{k} \sum_{l_j=0}^{q_j} \beta_{j,l_j} \Delta x_{j,t-l_j} + \lambda (y_{t-1} - \theta_0 - \sum_{j=1}^{k} \theta_j x_{j,t-1}) + \epsilon_t \]

Burada: - \(\Delta y_t\) ve \(\Delta x_{j,t}\) fark operatörü ile elde edilen değişkenlerin kısa dönem dinamiklerini gösterir. - \((y_{t-1} - \theta_0 - \sum_{j=1}^{k} \theta_j x_{j,t-1})\) uzun dönem denge ilişkisini ifade eder. - \(\lambda\) katsayısı, hata düzeltme teriminin katsayısıdır ve uzun dönem dengesine dönüş hızını gösterir.

Pesaran, Shin ve Smith (2001) tarafından geliştirilen sınır testi yöntemi, bu modelde \(\lambda\) katsayısının istatistiksel olarak anlamlı olup olmadığını test eder. Eğer \(\lambda\) anlamlı ise, modelde uzun dönem eşbütünleşme ilişkisi olduğu sonucuna varılabilir.

Son gösterim muhtemelen en ilginç olanı ve uygulamalı çalışmalarda genellikle en çok dikkat çekenidir. Buradaki amaç, tipik bir VAR çerçevesini koşullu hata düzeltme (CEC) formuna indirgemek suretiyle eşbütünleşmeyi test etmektir. Gerçekte, ilgilenilen CEC modeli, (3) numaralı denklemdeki modelle birebir karşılık gelen bir ARDL modelidir. Bunu görmek için, aşağıdaki 2. satırda \(y_t\) yerine denklem (1)’in sağ tarafını yerleştirin. Özellikle:

\[\Delta y_t = y_t - y_{t-1}\]

\[= a_0 + a_1 t + \psi^*(L)y_t + \sum_{j=1}^{k} \beta_j(L)x_{j,t} + \epsilon_t - y_{t-1}\]

\[= a_0 + a_1 t - y_{t-1} + \psi^*(1)y_{t-1} + \widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \beta_j(L)x_{j,t} + \epsilon_t\]

\[= a_0 + a_1 t - (1 - \psi^*(1)) y_{t-1} + \widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \beta_j(L) (x_{j,t-1} + \Delta x_{j,t}) + \epsilon_t\]

\[= a_0 + a_1 t - \psi(1)y_{t-1} + \widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \left( \beta_j(1) + (1 - L)\tilde{\beta}_j(L) \right) x_{j,t-1} + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[= a_0 + a_1 t - \psi(1)y_{t-1} + \sum_{j=1}^{k} \beta_j(1)x_{j,t-1}\]

\[+ \left(\widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \tilde{\beta}_j(L)\Delta x_{j,t-1} \right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t \tag{6}\]

Denklem (6), denklem (1)’deki ARDL modelinden türetilen CEC formudur. Bu denklemi şu şekilde yeniden yazabiliriz:

\[\Delta y_t = a_0 + a_1 t - \psi(1) \left( y_{t-1} - \sum_{j=1}^{k} \frac{\beta_j(1)}{\psi(1)} x_{j,t-1} \right)\]

\[+ \left(\widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \tilde{\beta}_j(L)\Delta x_{j,t-1} \right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[= a_0 + a_1 t - \psi(1)EC_{t-1}\]

\[+ \left(\widetilde{\psi^*}(L)\Delta y_{t-1} + \sum_{j=1}^{k} \tilde{\beta}_j(L)\Delta x_{j,t-1} \right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t \tag{7}\]

\(EC_t\) olarak gösterilen hata düzeltme teriminin, \(y_t\) ve \(x_{1,t}, \ldots, x_{k,t}\) eşbütünleşik olduğunda, aynı zamanda eşbütünleşme ilişkisi olduğu kolaylıkla doğrulanabilir. Gerçekte, PSS(2001) denklem (6)’nın (farklı gecikme değerlerinden soyutlanarak) VAR(p) modelinin CEC’i olduğunu göstermektedir:

\[\mathbf{\Phi}(L)(z_t - \boldsymbol{\mu} - \boldsymbol{\gamma}t) = \epsilon_t\]

Burada \(z_t\), \((y_t, x_{1,t}, \ldots, x_{k,t})^\top\) şeklinde bir \((k+1)\)-vektörüdür ve \(\boldsymbol{\mu}\) ve \(\boldsymbol{\gamma}\) sırasıyla sabit terim ve trend katsayılarının \((k+1)\)-vektörleridir, ve \(\mathbf{\Phi}(L) = I_{k+1} - \sum_{i=1}^{p} \mathbf{\Phi}_i L^i\) ise \((k+1)\) kare matris gecikme polinomudur. Bu özellikle önemlidir, çünkü CEC genellikle eşbütünleşmenin varlığını test etmek için bir platform olarak kullanılır.

Geleneksel olarak, Engle-Granger (1987), Phillips ve Ouliaris (1990) veya Johansen (1995) tarafından geliştirilen eşbütünleşme testleri, tipik olarak VAR modelindeki tüm değişkenlerin I(1) olmasını gerektirir. Bu açıkça, incelenen her bir değişkende birim kök varlığını test etmek için bir dizi ön test yapılmasını gerektirir ve özellikle birim kök testlerinin ilgilenilen birçok durumda büyüklük ve güç sorunlarından muzdarip olduğu bilindiğinden, yanlış sınıflandırmaya açıktır.

Buna karşılık, PSS(2001), ilgilenilen değişkenlerin I(0), I(1) veya karşılıklı olarak eşbütünleşik olup olmadığına karşı dayanıklı olmakla kalmayan, aynı zamanda yalnızca bilindik en küçük kareler regresyonlarında kullanılan tahmin ve çıkarımsal prosedürleri gerektirdiği için uygulaması önemli ölçüde daha kolay bir eşbütünleşme testi önermektedir.

Bu bağlamda, PSS(2001) ünlü sınır testini, CEC modelinin eşbütünleşme ilişkisindeki parametre anlamlılığına yönelik bir test olarak tartışmaktadır (6). Başka bir deyişle, test, aşağıdaki sıfır ve alternatif hipotezler için standart bir \(F\) veya Wald testidir:

\[H_0 : \psi(1) \cap \{\beta_j(1)\}_{j=1}^k = 0 \quad \text{(değişkenler eşbütünleşik değildir)}\]

\[H_A : \psi(1) \cup \{\beta_j(1)\}_{j=1}^k \neq 0 \quad \text{(değişkenler eşbütünleşiktir)}\]

Test istatistiği hesaplandıktan sonra, tüm değişkenlerin tamamen I(0) veya tamamen I(1) olduğu uç durumlara karşılık gelen iki asimptotik kritik değerle karşılaştırılır. Bu nedenle, bu kritik değerler, Brownian hareketlerinin integral fonksiyonlarını içeren standart olmayan bir karışım dağılımının sırasıyla alt ve üst kuyrukları içinde yer alır.

Test istatistiği alt kritik değerin altında olduğunda, sıfır hipotezi reddedilemez ve eşbütünleşmenin mümkün olmadığı sonucuna varılır. Buna karşılık, test istatistiği üst kritik değerin üzerinde olduğunda, sıfır hipotezi reddedilir ve eşbütünleşmenin gerçekten mümkün olduğu sonucuna varılır. Bu iki durumdan herhangi birinde, eşbütünleşme rankı hakkında bilgi sahibi olmak gerekli değildir.

Alternatif olarak, test istatistiği alt ve üst kritik değerler arasında kalırsa, test sonuçsuz kalır ve daha ileri gitmek için eşbütünleşme rankı hakkında bilgi gereklidir.

Burada ayrıca belirtmek gerekir ki, Narayan (2005)’te PSS (2001)’de sunulan asimptotik kritik değerlerin, \(T = 1000\) örnek büyüklükleri için türetildiklerinden dolayı çoğu pratik uygulama için gerçekçi olmadığı savunulmuştur. Buna göre, Narayan (2005), \(T = 30\)’dan \(T = 80\)’e kadar 5’er artışla değişen örnek büyüklükleri için kritik değerler sunmaktadır. Bunların, çoğu sonlu örnek ortamında çıkarımsal güvenilirliği iyileştireceğini savunmaktadırlar. EViews 10, kullanıcıya PSS (2001) veya Narayan (2005) kritik değerlerini kullanma seçeneği sunacaktır.

Burada ayrıca PSS(2001)’in CEC modelinin (6) beş alternatif yorumunu sunduğunu vurgulamak önemlidir. Bu yorumlar, deterministik terimlerin hata düzeltme terimine entegre olup olmadığına göre ayrılmaktadır. Deterministik terimler hata düzeltme terimine katkıda bulunduğunda, eşbütünleşme vektörünün açılımına dolaylı olarak yansıtılırlar. Başka bir deyişle, VAR(\(p\)) modelindeki \(\boldsymbol{\mu}\) ve \(\boldsymbol{\gamma}\) kısıtlıdır - eşbütünleşme vektöründeki elemanların doğrusal bir kombinasyonudur. Bu, açıkça, denklem (6)’daki \(a_0\) ve \(a_1\)’in de benzer şekilde kısıtlanması gerektiğini gösterir. Aşağıda, beş yorumun her biri için teorik (DGP) ve pratik regresyon (REG) modellerinin özetleri, uygun eşbütünleşme ilişkisi \((EC_t)\) ve sınır testi sıfır hipotezi \((H_0)\) ile birlikte verilmiştir.

Case 1: (Sabit Terim ve Trend Yok): \(a_0 = a_1 = 0\), yani, \((\boldsymbol{\mu} = \boldsymbol{\gamma} = \mathbf{0})\)

DGP:

\[\Delta y_t = -\psi(1)y_{t-1} + \sum_{j=1}^k \beta_j(1)x_{j,t-1}\]

\[+ \left(\psi^*(L)\Delta y_{t-1} + \sum_{j=1}^k \beta_j^*(L)\Delta x_{j,t-1}\right) + \sum_{j=1}^k \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[EC_t = y_t - \sum_{j=1}^k \frac{\beta_j(1)}{\psi(1)}x_{j,t}\]

\[H_0 : \psi(1) = \beta_j(1) = 0, \quad \forall j\]

Regresyon:

\[\Delta y_t = b_0 y_{t-1} + \sum_{j=1}^k b_j x_{j,t-1}\]

\[+ \sum_{i=1}^{p-1} c_{0,i}\Delta y_{t-i} + \sum_{j=1}^k \sum_{l_j=1}^{q_j-1} c_{j,l_j}\Delta x_{j,t-l_j} + \sum_{j=1}^k d_j \Delta x_{j,t} + \epsilon_t \quad (8)\]

\[EC_t = y_t - \sum_{j=1}^k \frac{b_j}{b_0}x_{j,t}\]

\[H_0 : b_0 = b_j = 0, \quad \forall j\]

Case 2: (Kısıtlı Sabit Terim ve Trend Yok): \(a_0 = \psi(1)\mu_y + \sum_{j=1}^k \beta_j(1)\mu_{x_j}\) ve \(a_1 = 0\) öyle ki \(\boldsymbol{\gamma} = \mathbf{0}\)

DGP:

\[\Delta y_t = -\psi(1) (y_{t-1} - \mu_y) + \sum_{j=1}^k \beta_j(1) (x_{j,t-1} - \mu_{x_j})\]

\[+ \left(\psi^*(L)\Delta y_{t-1} + \sum_{j=1}^k \beta_j^*(L)\Delta x_{j,t-1}\right) + \sum_{j=1}^k \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[EC_t = y_t - \mu_y - \sum_{j=1}^k \frac{\beta_j(1)}{\psi(1)}(x_{j,t} - \mu_{x_j})\]

\[H_0 : \psi(1) = \beta_j(1) = 0, \quad \forall j\]

Regresyon:

\[\Delta y_t = a_0 + b_0 y_{t-1} + \sum_{j=1}^k b_j x_{j,t-1}\]

\[+ \sum_{i=1}^{p-1} c_{0,i}\Delta y_{t-i} + \sum_{j=1}^k \sum_{l_j=1}^{q_j-1} c_{j,l_j}\Delta x_{j,t-l_j} + \sum_{j=1}^k d_j \Delta x_{j,t} + \epsilon_t \quad (9)\]

\[EC_t = y_t - \sum_{j=1}^k \frac{b_j}{b_0}x_{j,t} - \frac{a_0}{b_0}\]

\[H_0 : a_0 = b_0 = b_j = 0, \quad \forall j\]

Durum 3: (Kısıtlanmamış Sabit ve Trend Yok): \(a_0 \neq 0\) ve \(a_1 = 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(\boldsymbol{\gamma} = \mathbf{0}\)

DGP:

\[\Delta y_t = a_0 - \psi(1)y_{t-1} + \sum_{j=1}^{k} \beta_j(1)x_{j,t-1}\]

\[+ \left(\widetilde{\psi}^*(L)\Delta y_{t-1} + \sum_{j=1}^{k} \widetilde{\beta}_j(L)\Delta x_{j,t-1}\right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[EC_t = y_t - \sum_{j=1}^{k} \frac{\beta_j(1)}{\psi(1)} x_{j,t}\]

\[H_0 : \quad \psi(1) = \beta_j(1) = 0, \quad \forall j\]

Regresyon:

\[\Delta y_t = a_0 + b_0 y_{t-1} + \sum_{j=1}^{k} b_j x_{j,t-1}\]

\[+ \sum_{i=1}^{p-1} c_{0,i}\Delta y_{t-i} + \sum_{j=1}^{k} \sum_{l_j=1}^{q_j-1} c_{j,l_j}\Delta x_{j,t-l_j} + \sum_{j=1}^{k} d_j \Delta x_{j,t} + \epsilon_t \quad (10)\]

\[EC_t = y_t - \sum_{j=1}^{k} \frac{b_j}{b_0} x_{j,t}\]

\[H_0 : \quad b_0 = b_j = 0, \quad \forall j\]

Durum 4: (Kısıtlanmamış Sabit ve Kısıtlı Trend): \(a_0 \neq 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(a_1 = \psi(1)\gamma_y + \sum_{j=1}^{k} \beta_j(1)\gamma_{x_j}\)

DGP:

\[\Delta y_t = a_0 - \psi(1) \left(y_{t-1} - \gamma_y t\right) + \sum_{j=1}^{k} \beta_j(1) \left(x_{j,t-1} - \gamma_{x_j}t\right)\]

\[+ \left(\widetilde{\psi}^*(L)\Delta y_{t-1} + \sum_{j=1}^{k} \widetilde{\beta}_j(L)\Delta x_{j,t-1}\right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[EC_t = y_t - \gamma_y t - \sum_{j=1}^{k} \frac{\beta_j(1)}{\psi(1)} \left(x_{j,t} - \gamma_{x_j}t\right)\]

\[H_0 : \quad \psi(1) = \beta_j(1) = 0, \quad \forall j\]

Regresyon:

\[\Delta y_t = a_0 + a_1 t + b_0 y_{t-1} + \sum_{j=1}^{k} b_j x_{j,t-1}\]

\[+ \sum_{i=1}^{p-1} c_{0,i}\Delta y_{t-i} + \sum_{j=1}^{k} \sum_{l_j=1}^{q_j-1} c_{j,l_j}\Delta x_{j,t-l_j} + \sum_{j=1}^{k} d_j \Delta x_{j,t} + \epsilon_t \quad (11)\]

\[EC_t = y_t - \sum_{j=1}^{k} \frac{b_j}{b_0} x_{j,t} - \frac{a_1}{b_0} t\]

\[H_0 : \quad a_1 = b_0 = b_j = 0, \quad \forall j\]

Durum 5: (Kısıtlanmamış Sabit ve Kısıtlanmamış Trend): \(a_0 \neq 0\) ve \(a_1 \neq 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(\boldsymbol{\gamma} \neq \mathbf{0}\)

DGP:

\[\Delta y_t = a_0 + a_1 t - \psi(1)y_{t-1} + \sum_{j=1}^{k} \beta_j(1)x_{j,t-1}\]

\[+ \left(\widetilde{\psi}^*(L)\Delta y_{t-1} + \sum_{j=1}^{k} \widetilde{\beta}_j(L)\Delta x_{j,t-1}\right) + \sum_{j=1}^{k} \beta_j(L)\Delta x_{j,t} + \epsilon_t\]

\[EC_t = y_t - \sum_{j=1}^{k} \frac{\beta_j(1)}{\psi(1)} x_{j,t}\]

\[H_0 : \quad \psi(1) = \beta_j(1) = 0, \quad \forall j\]

Regresyon:

\[\Delta y_t = a_0 + a_1 t + b_0 y_{t-1} + \sum_{j=1}^{k} b_j x_{j,t-1}\]

\[+ \sum_{i=1}^{p-1} c_{0,i}\Delta y_{t-i} + \sum_{j=1}^{k} \sum_{l_j=1}^{q_j-1} c_{j,l_j}\Delta x_{j,t-l_j} + \sum_{j=1}^{k} d_j \Delta x_{j,t} + \epsilon_t \quad (12)\]

\[EC_t = y_t - \sum_{j=1}^{k} \frac{b_j}{b_0} x_{j,t}\]

\[H_0 : \quad b_0 = b_j = 0, \quad \forall j\]

Bu üç temsil, ARDL modellerinin zaman içindeki dinamiklerini, uzun dönem dengesini ve hata düzeltme mekanizmasını anlamak için kritik öneme sahiptir.

ARDL Model Uygulama Prosedürü

Burada yukarıdaki bölümleri pratikte nasıl uygulayacağımıza dair adım adım bir prosedür sunuyoruz.

1. Verinin doğasını anlamaya çalışın.
1. Tüm değişkenlerin I(d) mertebesinden entegre olduğundan ve d<2 olduğundan emin olun.
1. Deterministiklerin ARDL modeline nasıl gireceğini belirleyin. Bölüm 1 de özetlenen DGP 1,…,5 arasından seçim yapın.
1. Adım 3’te seçilen modelin uygun gecikme yapısını belirleyin.
1. Adım 4’teki modeli Olağan En Küçük Kareler (OLS) kullanarak tahmin edin.
1. Adım 5’ten elde edilen artıkların seri olarak ilişkisiz ve homoskedastik olduğundan emin olun.
1. Sınır Testini gerçekleştirin.
1. Uygun olduğunda ayarlama hızını tahmin edin.

Aşağıdaki akış şeması prosedürü göstermektedir.

Faiz Oranı Termin Yapısı (TSIR) ve Eşbütünleşme

Bu girdinin motivasyonu, klasik faiz oranı termin yapısı (TSIR) literatürüdür. Özetle, TSIR, farklı vadelere sahip tahvillerin getirileri arasında bir ilişki olduğunu varsayar.

Formel olarak:

\[ R(k,t) = \frac{1}{k} \sum_{j=1}^{k} E_t R(1,t+j-1) + L(k,t) \]

Burada \(E_t\), \(t\) zamanındaki bilgilere koşullu beklenti operatörüdür, \(R(k,t)\), \(t\) zamanında \(k\) dönemlik bir iskonto tahvilinin vadeye kadar getirisidir ve \(L(k,t)\) genellikle riski hesaba katan primlerdir. Eşbütünleşmenin mümkün olduğunu görmek için, aşağıdaki numara tekrar tekrar uygulanır:

\[ R(k,t) = R(k,t-1) + \Delta R(k,t), \]

Burada

\[ \Delta R(k,t) = R(k,t) - R(k,t-1), \]

Bu, aşağıdaki ifadeye yol açar:

\[ R(k,t) - R(1,t) = \frac{1}{k} \sum_{i=1}^{k-1} \sum_{j=1}^{i} E_t \Delta R(1,t+j) + L(k,t) \]

Şimdi, eğer \(R(k,t)\) I(1) süreçleriyse, \(\Delta R(1,t+j)\) I(0) süreçleri olmalıdır ve dolayısıyla \(L(k,t)\) de öyleyse, \(R(k,t) - R(1,t)\) doğrusal kombinasyonu I(0) süreçleridir. Başka bir deyişle, \(k\) dönemlik vadeye kadar getiri, her zaman birinci dönem vadeye kadar getiri ile eşbütünleşiktir ve eşbütünleşme vektörü \((1,-1)^T\)’dir. Aslında, biraz daha çalışma, bu ilkenin herhangi iki keyfi zaman \(k_1\) ve \(k_2\) arasındaki spread için de geçerli olduğunu gösterir. Yani,

\[ R(k_2,t) - R(k_1,t) = R(k_2,t) - R(1,t) + R(1,t) - R(k_1,t) \]

\[ = \frac{1}{k_2} \sum_{i=1}^{k_2-1} \sum_{j=1}^{i} E_t \Delta R(1,t+j) + L(k_2,t) - \frac{1}{k_1} \sum_{i=1}^{k_1-1} \sum_{j=1}^{i} E_t \Delta R(1,t+j) + L(k_1,t) \]

\[ \sim I(0) \] ### Uygulamaya Geçiş: Gerçek Verilerle Çalışma

Teorik temeli oluşturduğumuza göre, şimdi gerçek verilerle pratiğe geçiyoruz. Kanada olgunluk verileriyle çalışacağız. Bu veriler, Kanada Sosyoekonomik Veritabanı (Canadian Socioeconomic Database) veya kısaca CANSIM’den doğrudan toplanmıştır. Özellikle, iki tür pazarlanabilir borç enstrümanı arasındaki eşbütünleşme ilişkilerine odaklanacağız: Hazine Bonosu getirisi ve Benchmark Tahviller (diğer adıyla Hazine Notları) getirisi.

Hazine Bonosu Getirisi: Kısa vadeli (ihraç tarihinden itibaren 1 ay, 3 ay, 6 ay ve 1 yıl vadeli) iskonto güvencesi olan bir enstrümandır.
Benchmark Tahviller Getirisi: Orta vadeli (ihraç tarihinden itibaren 2 yıl, 5 yıl, 7 yıl ve 10 yıl vadeli) ve iki yılda bir faiz ödemesi yapan güvencelerdir.

Bu iki enstrüman arasındaki eşbütünleşme ilişkilerini inceleyerek, piyasa dinamiklerini daha iyi anlamayı hedefliyoruz.

library(readr)
gg <- read_csv("gg.csv")

# Gerekli kütüphaneleri yükle
library(ggplot2)
library(dplyr)
library(tidyr)

# Veri setini uzun formata dönüştür (DUM0708 hariç)
gg_long <- gg %>%
  pivot_longer(cols = c(BBY10Y, BBY2Y, BBY5Y, TBILL1M, TBILL1Y, TBILL3M, TBILL6M), 
               names_to = "Variable", values_to = "Value") %>%
  mutate(obs = as.Date(paste0(substr(obs, 1, 7), ":01"), format = "%Y:%m:%d"))  # Tarih formatını düzelt

# 2007 yılı için dikey çizgi eklemek üzere bir veri çerçevesi oluştur
vline_data <- data.frame(x = as.Date("2007-01-01"))

# Her bir değişken için ayrı grafikler oluştur
for (variable in unique(gg_long$Variable)) {
  p <- ggplot(gg_long %>% filter(Variable == variable), aes(x = obs, y = Value)) +
    geom_line() +
    geom_vline(data = vline_data, aes(xintercept = x), color = "red", linetype = "dashed") +
    labs(title = paste("Zaman Serisi Grafiği -", variable),
         x = "Tarih",
         y = "Değer") +
    theme_minimal()
  print(p)
}

Her grafik, Haziran 2007 civarında yapısal bir değişiklik sergiliyor ve bu, ABD konut krizinin başlangıcını işaret ediyor. Bu durumu, kırmızı dikey bir çizgi ile belirttik. Bu bilgiyi analizimize dahil edeceğiz ve kriz sonrası dönemi dum0708 adlı kukla değişkenle göstereceğiz. Bu değişken, Haziran 2007’yi takip eden her ay için 1 değerini alacak. Ayrıca, Kanada Merkez Bankası (CBC) üzerine yapılan küçük bir arka plan araştırması, Ocak 2001’den itibaren CBC’nin, 90’ların sonundaki dot-com çöküşünden ve o on yılın başlarındaki deflasyonist dönemden kurtulmak için yeni bir şeffaflık ve enflasyon hedefleme önlemleri setine bağlı kalacağını ortaya koyuyor. Bu nedenle, çok fazla politika paradigma değişikliğini analiz etmek zorunda kalmamak için yalnızca Ocak 2001 sonrası dönemdeki verilere odaklanacağız.

Analizimize, incelenen serilerden hiçbirinin 2. dereceden veya daha yüksek dereceden entegre olmadığından emin olarak başlıyoruz. Bunu yapmak için, her bir serinin birinci farkı üzerinde bir birim kök testi uyguluyoruz. Bu durumda, standart ADF (Augmented Dickey-Fuller) testi yeterli olacaktır.

Bu test, serilerin durağanlık özelliklerini belirlemek için kullanılır. Eğer bir serinin birinci farkı durağan çıkarsa, bu seri I(1) olarak kabul edilir ve daha yüksek dereceden entegre olmadığı anlaşılır. Bu adım, eşbütünleşme analizi için gerekli bir ön koşuldur.

# Gerekli kütüphaneleri yükle
library(tseries)

# Serilerin birinci farklarını al ve ADF testi uygula
adf_results <- list()

# BBY10Y için ADF testi
adf_results$BBY10Y <- adf.test(diff(gg$BBY10Y))

# BBY2Y için ADF testi
adf_results$BBY2Y <- adf.test(diff(gg$BBY2Y))

# BBY5Y için ADF testi
adf_results$BBY5Y <- adf.test(diff(gg$BBY5Y))

# TBILL1M için ADF testi
adf_results$TBILL1M <- adf.test(diff(gg$TBILL1M))

# TBILL1Y için ADF testi
adf_results$TBILL1Y <- adf.test(diff(gg$TBILL1Y))

# TBILL3M için ADF testi
adf_results$TBILL3M <- adf.test(diff(gg$TBILL3M))

# TBILL6M için ADF testi
adf_results$TBILL6M <- adf.test(diff(gg$TBILL6M))

# Test sonuçlarını görüntüle
adf_results

## $BBY10Y
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$BBY10Y)
## Dickey-Fuller = -7.3075, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $BBY2Y
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$BBY2Y)
## Dickey-Fuller = -7.0347, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $BBY5Y
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$BBY5Y)
## Dickey-Fuller = -7.6731, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $TBILL1M
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$TBILL1M)
## Dickey-Fuller = -4.1364, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $TBILL1Y
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$TBILL1Y)
## Dickey-Fuller = -5.8716, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $TBILL3M
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$TBILL3M)
## Dickey-Fuller = -4.3605, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $TBILL6M
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(gg$TBILL6M)
## Dickey-Fuller = -5.114, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Her bir seri için p-değeri 0 olduğundan ve boş hipotez (null hypothesis) bir birim kök (unit root) olduğundan, tüm anlamlılık düzeylerinde boş hipotezi reddediyoruz. Özellikle, test birinci farklar (first differences) üzerinde uygulandığı için, birinci farklarda birim kök olmadığı sonucuna varıyoruz. Bu nedenle, her bir seri ya I(0) (durağan) ya da I(1) (birinci dereceden entegre) olmalıdır.

Bu sonuç, analizimizin ikinci adımına geçebileceğimizi gösteriyor. Bir sonraki adım, seriler arasında eşbütünleşme (cointegration) ilişkisi olup olmadığını incelemek olacaktır. Bu, özellikle serilerin uzun dönemli bir denge ilişkisi içinde olup olmadığını belirlemek için önemlidir.

Deterministik Belirlemeler

Verilere uygun bir model seçmek hem sanat hem de bilimdir. Bununla birlikte, birkaç rehber ilke bulunmaktadır. Serilerin sıfır etrafında merkezlenmediği herhangi bir model genellikle bir sabit terim gerektirirken, serilerde bir eğilim (trend) bulunan modeller genellikle bir eğilim terimi içerdiğinde daha iyi uyum sağlar. Bu dersin 1. ve 2. bölümlerinde, Case 1’den Case 5’e kadar beş farklı Veri Üretim Süreci (DGP) belirlemesinden seçim yapma olasılığı tartışıldı. Aslında, çeşitli değişken kombinasyonlarıyla birkaç farklı model belirlemesini ele alacağız.

Model 1:

Amaç: 10 Yıllık Benchmark Tahvil Getirisi (10 Year Benchmark Bond Yield) ile 1 Aylık Hazine Bonosu Getirisi (1 Month T-Bill) arasındaki ilişkiyi incelemek.

Sabit Terim: Sabit terim, eşbütünleşme ilişkisine dahil edilecek şekilde kısıtlanmıştır.
DGP ve Regresyon Modeli: Bu model, Case 2’de belirtilen DGP ve regresyon modeline karşılık gelir.
Model 2:

Amaç: 6 Aylık, 3 Aylık ve 1 Aylık Hazine Bonoları (6, 3, and 1 Month T-Bills) arasındaki ilişkiyi incelemek.

Sabit Terim: Sabit terim serbest bırakılmıştır (unrestricted).
DGP ve Regresyon Modeli: Bu model, Case 3’te belirtilen DGP ve regresyon modeline karşılık gelir.
Model 3:

Amaç: 2 Yıllık Benchmark Tahvil Getirisi (2 Year Benchmark Bond Yield) ile 1 Yıllık ve 1 Aylık Hazine Bonoları (1 Year and 1 Month T-Bills) arasındaki ilişkiyi incelemek.

Sabit Terim: Sabit terim serbest bırakılmıştır (unrestricted).
DGP ve Regresyon Modeli: Bu model, Case 3’te belirtilen DGP ve regresyon modeline karşılık gelir.

ARDL Gecikme Yapısının Belirlenmesi

Model için uygun gecikme sayısını seçmek, yine hem bilim hem de sanattır. Ekonomik teori tarafından gecikme sayısı belirtilmediği sürece, ekonometristin optimal gecikme uzunluğunu seçmek için çeşitli araçları vardır. Bir olasılık, bağımlı değişken için maksimum gecikme sayısını (örneğin, p) ve her bir açıklayıcı değişken için maksimum gecikme sayısını (örneğin, q belirlemek ve ardından bu belirlemeyi kullanarak oluşturulabilecek tüm farklı gecikme kombinasyonlarıyla bir dizi regresyon çalıştırmaktır. Özellikle, k açıklayıcı değişken varsa, {1,…,p} kümesi ve k ek {0,…,q} kümesi kullanılarak oluşturulabilecek maksimum kombinasyon sayısı \(p×(q+1)^k\) olacaktır. Örneğin, varsayılan değerleri olan p=q=4 ile, dikkate alınacak toplam model sayısı 100 olacaktır. Optimal kombinasyon, daha sonra Akaike (AIC), Schwarz (BIC), Hannan-Quinn (HQ) veya hatta düzeltilmiş R kare gibi bir bilgi kriterini en aza indiren kombinasyon olarak belirlenir.

Model 1: Eşbütünleşme yok (No Cointegrating Relationship)

library(ARDL)

model <- auto_ardl(BBY10Y ~ TBILL1M | DUM0708, data=gg, max_order = 4)

model$best_order

##  BBY10Y TBILL1M 
##       1       3

bm <- model$best_model
summary(bm)

## 
## Time series regression with "ts" data:
## Start = 4, End = 196
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37474 -0.12581 -0.02084  0.12194  0.52703 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.18349    0.08674   2.115 0.035724 *  
## L(BBY10Y, 1)   0.95167    0.01907  49.901  < 2e-16 ***
## TBILL1M        0.25229    0.06991   3.609 0.000395 ***
## L(TBILL1M, 1) -0.38160    0.09398  -4.060 7.21e-05 ***
## L(TBILL1M, 2)  0.01905    0.09393   0.203 0.839539    
## L(TBILL1M, 3)  0.12068    0.06739   1.791 0.074948 .  
## DUM0708       -0.09875    0.05089  -1.940 0.053836 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1888 on 186 degrees of freedom
## Multiple R-squared:  0.9803, Adjusted R-squared:  0.9796 
## F-statistic:  1539 on 6 and 186 DF,  p-value: < 2.2e-16

library(lmtest)
library(strucchange)

bgtest(bm, order=2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  bm
## LM test = 0.8031, df = 2, p-value = 0.6693

Boş hipotez (null hypothesis), artıkların serisel olarak ilişkisiz (serially uncorrelated) olduğu yönündedir. F-istatistiği için elde edilen p-değeri 0.6693 olduğundan, bu boş hipotezi reddedemeyiz. Bu nedenle, artıkların serisel olarak ilişkisiz olduğu sonucuna varırız.

bptest(bm)

## 
##  studentized Breusch-Pagan test
## 
## data:  bm
## BP = 10.297, df = 6, p-value = 0.1127

Boş hipotez (null hypothesis), artıkların homoskedastik (homoskedastic) olduğu yönündedir. F-istatistiği için elde edilen p-değeri 0.1127 olduğundan, bu boş hipotezi %10 anlamlılık düzeyinde bile reddedemeyiz. Bu nedenle, artıkların %10 anlamlılık düzeyinde homoskedastik olduğu sonucuna varırız.

boundtest <- bounds_f_test(bm, case = 2, alpha = 0.05)
boundtest$tab

##   statistic Lower-bound I(0) Upper-bound I(1) alpha  p.value
## F  2.275448         3.620887         4.135148  0.05 0.334749

boundtest$PSS2001parameters

## Lower-bound I(0) Upper-bound I(1) 
##             3.62             4.16

multipliers(bm)

##          Term  Estimate Std. Error   t value     Pr(>|t|)
## 1 (Intercept) 3.7970471  1.0621823 3.5747602 0.0004466781
## 2     TBILL1M 0.2156494  0.3290297 0.6554102 0.5130134921

F-istatistiği değeri 2.275448, açıkça I(0) kritik değer sınırının altındadır. Bu serinin 2. bölümündeki analizimiz, dengeleyici bir ilişki olmadığı yönündeki boş hipotezi reddedemediğimizi göstermektedir.

Aslında, uzun dönem denkleminin uyumunu ve bağımlı değişkeni görselleştirebiliriz. Bunu, EC (hata düzeltme) terimini çıkarıp bağımlı değişkenden çıkararak yapabiliriz. Bu işlem şu şekilde gerçekleştirilir:

ce <- coint_eq(bm, case = 2)

gg <- gg %>% mutate(obs = as.Date(paste0(substr(obs, 1, 7), ":01"), format = "%Y:%m:%d"))

ggplot(data = gg, aes(x = obs)) +
  geom_line(aes(y = BBY10Y, color = "BBY10Y")) +
  geom_line(aes(y = ce, color = "Cointegrating Eq.")) +
  labs(title = "BBY10Y vs Cointegrating Equation",
       x = "Date",
       y = "Value",
       color = "Series") +
  theme_minimal()

normalize_zscore <- function(x) {
  (x - mean(x)) / sd(x)
}

# BBY10Y ve ce serilerini normalize et
gg$BBY10Y_norm <- normalize_zscore(gg$BBY10Y)
gg$ce_norm <- normalize_zscore(ce)

ggplot(data = gg, aes(x = obs)) +
  geom_line(aes(y = BBY10Y_norm, color = "BBY10Y (Normalized)")) +
  geom_line(aes(y = ce_norm, color = "Cointegrating Eq. (Normalized)")) +
  labs(title = "Normalized BBY10Y vs Cointegrating Equation",
       x = "Date",
       y = "Normalized Value",
       color = "Series") +
  theme_minimal()

## Model 2: Usual Cointegrating Relationship

Bu modelde, bağımlı değişken 6 Aylık Hazine Bonosu (6 Months T-Bill)’dir. Dinamik açıklayıcı değişkenler ise 3 Aylık ve 1 Aylık Hazine Bonoları (3 and 1 Month T-Bills)’dir. Ayrıca, modelde kısıtlanmamış bir sabit terim (unrestricted constant) bulunmaktadır, bu da Case 3’e karşılık gelir. Modelde, dinamik olmayan açıklayıcı değişken olarak dum0708 (2008 finansal krizi sonrası dönemi temsil eden kukla değişken) kullanılmıştır.

model2 <- auto_ardl(TBILL6M ~ TBILL1M + TBILL3M | DUM0708, data=gg, max_order = 4)

model2$best_order

## TBILL6M TBILL1M TBILL3M 
##       4       4       4

model2bm <- model2$best_model

summary(model2bm)

## 
## Time series regression with "ts" data:
## Start = 5, End = 196
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.210868 -0.027454 -0.003247  0.024063  0.188327 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.004293   0.018920  -0.227  0.82076    
## L(TBILL6M, 1)  0.566397   0.079792   7.098 2.98e-11 ***
## L(TBILL6M, 2)  0.049489   0.091513   0.541  0.58934    
## L(TBILL6M, 3) -0.074471   0.095175  -0.782  0.43499    
## L(TBILL6M, 4) -0.160658   0.078261  -2.053  0.04157 *  
## TBILL1M       -0.417253   0.076276  -5.470 1.52e-07 ***
## L(TBILL1M, 1) -0.187524   0.082127  -2.283  0.02361 *  
## L(TBILL1M, 2)  0.236562   0.082243   2.876  0.00452 ** 
## L(TBILL1M, 3) -0.258513   0.082340  -3.140  0.00198 ** 
## L(TBILL1M, 4)  0.056752   0.076760   0.739  0.46069    
## TBILL3M        1.246926   0.076003  16.406  < 2e-16 ***
## L(TBILL3M, 1) -0.170158   0.116036  -1.466  0.14432    
## L(TBILL3M, 2) -0.312487   0.117784  -2.653  0.00871 ** 
## L(TBILL3M, 3)  0.292106   0.120473   2.425  0.01633 *  
## L(TBILL3M, 4)  0.135094   0.119804   1.128  0.26101    
## DUM0708        0.018712   0.015467   1.210  0.22798    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06144 on 176 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.9978 
## F-statistic:  5867 on 15 and 176 DF,  p-value: < 2.2e-16

bgtest(model2bm, order=2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  model2bm
## LM test = 12.986, df = 2, p-value = 0.001514

bptest(model2bm)

## 
##  studentized Breusch-Pagan test
## 
## data:  model2bm
## BP = 79.239, df = 15, p-value = 9.619e-11

Her iki testten elde edilen p-değerlerine göre, her iki testte de boş hipotezi reddedeceğiz. Açıkça görülüyor ki, hem serisel korelasyon (serial correlation) hem de heteroskedastisite (heteroskedasticity) sorunumuz var. İlk sorunu çözmek için, hem bağımlı değişken hem de açıklayıcı değişkenler için gecikme sayısını artıracağız. İkinci sorunu çözmek için ise, HAC (Newey-West) kovaryans matrisi düzeltmesi kullanacağız. Bu düzeltme, tahmin sırasında hesaplanan test istatistiklerinin değerlerini düzeltir. Unutmayın ki, serisel korelasyon yanlı (biased) sonuçlara yol açabilirken, heteroskedastisite sadece etkin olmayan (inefficient) tahminlere neden olur. Bu nedenle, serisel korelasyonu ortadan kaldırmak birincil öneme sahiptir. Bu iki görevi de şimdi gerçekleştireceğiz.

model2 <- auto_ardl(TBILL6M ~ TBILL1M + TBILL3M | DUM0708, data=gg, max_order = 8)

model2$best_order

## TBILL6M TBILL1M TBILL3M 
##       7       3       6

model2bm <- model2$best_model

summary(model2bm)

## 
## Time series regression with "ts" data:
## Start = 8, End = 196
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.146752 -0.030125 -0.001824  0.023595  0.197798 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.01569    0.01805  -0.869  0.38607    
## L(TBILL6M, 1)  0.53528    0.07783   6.878 1.13e-10 ***
## L(TBILL6M, 2)  0.06139    0.08636   0.711  0.47811    
## L(TBILL6M, 3) -0.05789    0.08915  -0.649  0.51701    
## L(TBILL6M, 4) -0.23298    0.08721  -2.672  0.00829 ** 
## L(TBILL6M, 5) -0.02242    0.08707  -0.257  0.79715    
## L(TBILL6M, 6)  0.09285    0.07869   1.180  0.23965    
## L(TBILL6M, 7)  0.08341    0.02808   2.970  0.00341 ** 
## TBILL1M       -0.38353    0.07316  -5.243 4.68e-07 ***
## L(TBILL1M, 1) -0.13295    0.07975  -1.667  0.09733 .  
## L(TBILL1M, 2)  0.21496    0.08058   2.668  0.00838 ** 
## L(TBILL1M, 3) -0.24937    0.07798  -3.198  0.00165 ** 
## TBILL3M        1.23554    0.07246  17.052  < 2e-16 ***
## L(TBILL3M, 1) -0.20122    0.11351  -1.773  0.07808 .  
## L(TBILL3M, 2) -0.28649    0.11518  -2.487  0.01384 *  
## L(TBILL3M, 3)  0.26050    0.11473   2.271  0.02444 *  
## L(TBILL3M, 4)  0.18692    0.06864   2.723  0.00715 ** 
## L(TBILL3M, 5)  0.07669    0.06998   1.096  0.27474    
## L(TBILL3M, 6) -0.17659    0.06147  -2.873  0.00459 ** 
## DUM0708        0.02116    0.01470   1.439  0.15194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05784 on 169 degrees of freedom
## Multiple R-squared:  0.9982, Adjusted R-squared:  0.998 
## F-statistic:  4948 on 19 and 169 DF,  p-value: < 2.2e-16

bgtest(model2bm, order=2)

## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  model2bm
## LM test = 0.036238, df = 2, p-value = 0.982

library(sandwich)
library(lmtest)

# Apply HAC (Newey-West) standard errors
bm_tolm <- to_lm(model2bm, fix_names = TRUE)
hac_se <- NeweyWest(bm_tolm)
coeftest(bm_tolm, vcov = hac_se)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.015689   0.016770 -0.9355 0.3508546    
## L.TBILL6M.1  0.535285   0.107160  4.9952 1.456e-06 ***
## L.TBILL6M.2  0.061393   0.083659  0.7338 0.4640597    
## L.TBILL6M.3 -0.057886   0.116375 -0.4974 0.6195472    
## L.TBILL6M.4 -0.232976   0.093070 -2.5032 0.0132557 *  
## L.TBILL6M.5 -0.022417   0.059555 -0.3764 0.7070911    
## L.TBILL6M.6  0.092852   0.051861  1.7904 0.0751798 .  
## L.TBILL6M.7  0.083409   0.024765  3.3680 0.0009378 ***
## TBILL1M     -0.383531   0.092330 -4.1539 5.179e-05 ***
## L.TBILL1M.1 -0.132955   0.048983 -2.7143 0.0073307 ** 
## L.TBILL1M.2  0.214963   0.106278  2.0227 0.0446853 *  
## L.TBILL1M.3 -0.249368   0.171770 -1.4518 0.1484230    
## TBILL3M      1.235544   0.111623 11.0689 < 2.2e-16 ***
## L.TBILL3M.1 -0.201225   0.052372 -3.8422 0.0001723 ***
## L.TBILL3M.2 -0.286489   0.080481 -3.5597 0.0004822 ***
## L.TBILL3M.3  0.260496   0.232192  1.1219 0.2634966    
## L.TBILL3M.4  0.186918   0.068976  2.7099 0.0074242 ** 
## L.TBILL3M.5  0.076687   0.049132  1.5608 0.1204328    
## L.TBILL3M.6 -0.176593   0.054349 -3.2492 0.0013962 ** 
## DUM0708      0.021155   0.012790  1.6541 0.0999674 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bptest(bm_tolm)

## 
##  studentized Breusch-Pagan test
## 
## data:  bm_tolm
## BP = 78.633, df = 19, p-value = 3.195e-09

F-istatistiği p-değeri 0.691, artık serisel korelasyon (serial correlation) sorunu olmadığını gösteriyor.

Eşbütünleşme (cointegration) varlığını test etmek için, tekrar Uzun Dönem Formu ve Sınır Testi (Long Run Form and Bounds Test) görünümüne geçiyoruz.

boundtest2 <- bounds_f_test(model2bm, case = 3, alpha = 0.05)
boundtest2$tab

##   statistic Lower-bound I(0) Upper-bound I(1) alpha      p.value
## F  9.018007         3.790645         4.806771  0.05 0.0007188644

bound_t_test2 <- bounds_t_test(model2bm, case = 3, alpha = 0.05)
bound_t_test2$tab

##   statistic Lower-bound I(0) Upper-bound I(1) alpha      p.value
## t -5.062143         -2.86012        -3.507306  0.05 0.0006556881

multipliers(model2bm)

##          Term    Estimate Std. Error    t value     Pr(>|t|)
## 1 (Intercept) -0.02903507 0.03259302 -0.8908372 3.742830e-01
## 2     TBILL1M -1.01952505 0.15737102 -6.4784805 9.726465e-10
## 3     TBILL3M  2.02712791 0.15561964 13.0261701 2.654778e-27

F-istatistiği değeri 16.94, açıkça I(1) kritik değer sınırından büyüktür. Bu serinin 2. bölümündeki analizimiz, dengeleyici bir ilişki olmadığı yönündeki boş hipotezi reddettiğimizi göstermektedir. Ayrıca, boş hipotezi reddettiğimiz ve eşbütünleşme ilişkisine bir sabit veya eğilim terimi eklemediğimiz için, bu serinin 2. bölümündeki açıklamalarımız, hangi alternatifin geçerli olduğunu belirlemek için t-Sınır Testi (t-Bounds Test) kritik değerlerini kullanabileceğimizi göstermektedir. Bu özel durumda, t-istatistiğinin mutlak değeri |-5.062143 | = 5.062143 ’dir ve bu değer, hem I(0) hem de I(1) t-sınırlarının mutlak değerinden büyüktür. Hatırlayacağınız gibi, bu durum, t-Sınır Testi boş hipotezini reddetmemiz ve eşbütünleşme ilişkisinin ya olağan türde olduğu ya da geçerli ancak dejenere olduğu sonucuna varmamız gerektiğini gösterir. Bununla birlikte, bağımlı değişken ile dengeleyici denklem arasındaki uyuma baktığımızda, ilişkinin gerçekten geçerli olduğuna inanmamız gerekir. Grafik aşağıda sunulmuştur.

ce2 <- coint_eq(model2bm, case = 3)

ggplot(data = gg, aes(x = obs)) +
  geom_line(aes(y = TBILL6M, color = "TBILL6M")) +
  geom_line(aes(y = ce2, color = "Cointegrating Eq.")) +
  labs(title = "TBILL6M vs Cointegrating Equation",
       x = "Date",
       y = "Value",
       color = "Series") +
  theme_minimal()

# TBILL6M ve ce serilerini normalize et
gg$TBILL6M_norm <- normalize_zscore(gg$TBILL6M)
gg$ce2_norm <- normalize_zscore(ce2)

ggplot(data = gg, aes(x = obs)) +
  geom_line(aes(y = TBILL6M_norm, color = "TBILL6M (Normalized)")) +
  geom_line(aes(y = ce2_norm, color = "Cointegrating Eq. (Normalized)")) +
  labs(title = "Normalized TBILL6M vs Cointegrating Equation",
       x = "Date",
       y = "Normalized Value",
       color = "Series") +
  theme_minimal()

Bu özel durumda, ayarlama hızı denklemini (speed of adjustment equation) incelemek mantıklıdır.

summary(recm(model2bm, case = 3))

## 
## Time series regression with "zooreg" data:
## Start = 8, End = 196
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.146752 -0.030125 -0.001824  0.023595  0.197798 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.015689   0.007515  -2.088  0.03831 *  
## d(L(TBILL6M, 1))  0.075625   0.099193   0.762  0.44687    
## d(L(TBILL6M, 2))  0.137018   0.092541   1.481  0.14055    
## d(L(TBILL6M, 3))  0.079132   0.083962   0.942  0.34728    
## d(L(TBILL6M, 4)) -0.153844   0.078439  -1.961  0.05147 .  
## d(L(TBILL6M, 5)) -0.176261   0.070487  -2.501  0.01334 *  
## d(L(TBILL6M, 6)) -0.083409   0.027581  -3.024  0.00288 ** 
## d(TBILL1M)       -0.383531   0.069571  -5.513 1.28e-07 ***
## d(L(TBILL1M, 1))  0.034405   0.094151   0.365  0.71525    
## d(L(TBILL1M, 2))  0.249368   0.075218   3.315  0.00112 ** 
## d(TBILL3M)        1.235544   0.069000  17.907  < 2e-16 ***
## d(L(TBILL3M, 1)) -0.061020   0.168417  -0.362  0.71757    
## d(L(TBILL3M, 2)) -0.347508   0.133906  -2.595  0.01028 *  
## d(L(TBILL3M, 3)) -0.087013   0.081725  -1.065  0.28851    
## d(L(TBILL3M, 4))  0.099906   0.073046   1.368  0.17320    
## d(L(TBILL3M, 5))  0.176593   0.060766   2.906  0.00414 ** 
## DUM0708           0.021155   0.010073   2.100  0.03719 *  
## ect              -0.540340   0.103275  -5.232 4.86e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0575 on 171 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.9138, Adjusted R-squared:  0.9053 
## F-statistic: 106.7 on 17 and 171 DF,  p-value: < 2.2e-16

Beklendiği gibi, EC terimi (hata düzeltme terimi), burada ect olarak temsil edilmektedir, negatiftir ve ilişkili katsayı tahmini -0.540340’dır. Bu, herhangi bir dengesizlik hareketinin yaklaşık %54.03’nün bir dönem içinde düzeltildiği anlamına gelir. Ayrıca, çok büyük bir t-istatistiği (yani -5.232) göz önüne alındığında, bu katsayının oldukça anlamlı olduğu sonucuna da varabiliriz. Daha fazla ayrıntı için bu serinin 2. bölümüne bakınız.

delay_15 <- multipliers(model2bm, type = 15, se=TRUE)
plot_delay(delay_15, interval = 0.95)

Gecikme çarpanları, etkinin zaman içinde nasıl değiştiğini gösterir.

Örneğin, etki başlangıçta daha güçlü olabilir ve zamanla azalabilir veya sabit kalabilir.

Pozitif çarpanların zaman içinde nasıl değiştiği, etkinin kalıcılığı veya geçiciliği hakkında bilgi verir.

TBILL3M’in pozitif gecikme çarpanları, 3 Aylık Hazine Bonosu getirisindeki artışların 6 Aylık Hazine Bonosu getirisini olumlu yönde etkilediğini gösterir.

Bu etki, zaman içinde azalıyor ve grafikteki güven aralıkları, etkinin istatistiksel olarak anlamlı olup olmadığını belirler.

Bu bulgu, kısa vadeli faiz hareketlerinin orta vadeli enstrümanlar üzerindeki etkisini anlamak için önemlidir.

bm_lm <- to_lm(model2bm, data_class = "ts")

resettest(bm_lm, power = 2)

## 
##  RESET test
## 
## data:  bm_lm
## RESET = 2.372, df1 = 1, df2 = 168, p-value = 0.1254

p-değeri 0.1254 > 0.05 olduğu için, boş hipotezi reddedemeyiz. Bu, modelin fonksiyonel formunun doğru belirlendiği anlamına gelir (yani, modelde eksik değişken veya yanlış fonksiyonel form yoktur).

jarque.bera.test(residuals(model2bm))

## 
##  Jarque Bera Test
## 
## data:  residuals(model2bm)
## X-squared = 29.425, df = 2, p-value = 4.079e-07

p-değeri 0.05’ten çok küçük olduğu için, boş hipotezi reddederiz. Bu, artıkların normal dağılmadığını gösterir. Normal dağılım varsayımının ihlali, özellikle hipotez testleri ve güven aralıkları için sorun yaratabilir.

fluctuations <- efp(bm_tolm$full_formula, data = bm_tolm$model)

sctest(fluctuations)

## 
##  Recursive CUSUM test
## 
## data:  fluctuations
## S = 0.86565, p-value = 0.08993

p-değeri 0.08993 > 0.05 olduğu için, boş hipotezi reddedemeyiz. Bu, modelin parametrelerinin zaman içinde sabit olduğunu gösterir (yani, yapısal kırılma yoktur). Ancak, p-değeri 0.05’e yakın olduğu için, parametrelerde hafif bir değişim olabileceğine dair bir şüphe olabilir.

plot(fluctuations)

NARDL: A Study of Bosnia’s Tourism Sector

Turizm, birçok ekonomi için hayati bir gelir kaynağıdır. Bosna-Hersek örneğinde ise, olağanüstü tarihi ve doğal turistik cazibesine rağmen, turizm 1990’ların başındaki korkunç saldırı döneminden önce ülkenin GSYH büyümesine kayda değer bir katkı sağlamıyordu. O dönemde Bosna ekonomisi büyük ölçüde doğal kaynakların (özellikle metal cevherleri ve ormancılık) kullanımı, hidroelektrik üretimi ve imalat sanayisine dayanıyordu. Bosna, 1996’da saldırıların sona ermesinin ardından bu endüstrilerin bir kısmını yeniden canlandırmayı başarsa da, sanayi üretimi gibi oldukça umut vadeden bazı sektörleri eski canlılığına kavuşturmakta zorlandı. Bu sektörlerin yerine turizm, Bosna’nın GSYH gelişimine giderek daha önemli bir katkı sağlayan bir alan haline gelmiştir.

Bosna’ya yönelik yabancı turist ilgisinin artışının büyük bir kısmı, aktif pazarlama kampanyalarından ve Avrupa Birliği’nin teşviklerinden kaynaklanmaktadır (bkz. AB’nin 2007/145-210 Projesi). 1984 Kış Olimpiyatları’na ev sahipliği yapmış olan Bosna’nın başkenti Saraybosna, aynı zamanda geniş çapta popüler olan Saraybosna Film Festivali’ne de ev sahipliği yapmaktadır. Güneyde yer alan ve UNESCO tarafından koruma altına alınmış Mostar kenti ise Red Bull Cliff Diving Dünya Serisi’ne ev sahipliği yapmakta olup, Avrupa’da ziyaret edilmesi gereken en güzel şehirler listelerinde sıkça yer almaktadır. Ayrıca, 1990’ların başındaki kanlı geçmişine rağmen ülke bugün, küresel çapta büyük ekonomilerde yaşayan ve çalışan 2,2 milyon kişilik (saldırı öncesi nüfusunun %55’i) bir diasporaya sahiptir. Bu faktörler, ülkenin görece düşük yaşam maliyeti ile birleştiğinde, Bosna-Hersek’i daha büyük ekonomilerden gelen turistler için oldukça cazip bir destinasyon haline getirmektedir.

Öte yandan, Bosna’nın iç turizm sektörü, yabancı turizme kıyasla geride kalmasına rağmen, benzer bir canlanma yaşamıştır. Ülkenin hem yaz hem de kış turizmine elverişli olması nedeniyle, Adriyatik Denizi kıyısındaki sahil şehirlerinde yapılan plaj aktiviteleri ile orta batı bölgesindeki kış kayak turizmi, yerel halk için yıl boyunca ülkelerini keşfetme fırsatı sunmaktadır. İç turizm, ayrıca sayısı yaklaşık 200’ü bulan yerel doğa yürüyüşü ve dağcılık kulüpleri tarafından teşvik edilmekte olup, bu kulüpler giderek daha fazla yabancı turistlere yönelik çok günlük rehberli turlar sunmaya başlamaktadır. Son olarak, Bosna’nın, saldırı dönemi boyunca ülkede kalmış önemli bir orta yaş grubu nüfusuna sahip olduğu da belirtilmelidir. Bu demografik grubun, savaş sonrası toparlanma yıllarında ekonomik fırsatları büyük ölçüde sınırlı kalmış olup, günümüzde uluslararası seyahatleri karşılayamayacak durumda olan birçok kişi, bunun yerine yurt içi gezilerle yetinmektedir.

Aslında, turizmin ekonomik büyüme üzerinde olumlu bir etkisi olduğunu gösteren çok sayıda çalışma bulunmaktadır. Bu araştırma alanını yönlendiren temel fikir, turizmin altyapı gelişimini, doğrudan yabancı yatırımları ve modern internet ve mobil bağlantı ağlarını teşvik eden olumlu bir dışsallık olduğudur. Buradaki amacımız, Bosna’nın turizm sektörünün genel ekonomik gelişimi üzerindeki büyüklüğünü ve asimetrilerini incelemektir.

Geleneksel olarak, bu dinamikler vektör hata düzeltme (VEC) modelleri ile incelenebilir. Bununla birlikte, bu model sınıfı, sistem değişkenlerinin en az birinci dereceden bütünleşik olduğunu varsayar ve birden fazla eşit derecede olası eşbütünleşme ilişkisini dışlamaz. Alternatif olarak, klasik ARDL modelleri, sistem değişkenleri arasında farklı bütünleşme derecelerine izin verir (maksimum bütünleşme derecesi 2’den küçük olmak kaydıyla) ve değişkenler arasında tekil ve benzersiz bir eşbütünleşme ilişkisi olduğunu varsayar. Ancak, bu çerçeve, uzun vadeli (eşbütünleşik) ilişkinin regresörlerin simetrik lineer bir kombinasyonu olduğunu varsaymaktadır.

Klasik ARDL çerçevesi birçok uygulama için oldukça makul olsa da, davranışsal finans ve ekonomi literatüründeki doğrusal olmayanlık ve asimetri konularını ele alamaz; ki bu durumlar da pratikte sıkça karşılaşılmaktadır (bkz. Kahneman ve Tversky, 1979; Shiller, 2005). Bu sınırlamayı aşmak için Shin, Yu ve Greenwood-Nimmo (2014), kısa ve uzun vadeli doğrusal olmayanlıkları, dağıtılmış gecikmeli değişkenlerin pozitif ve negatif kısmi toplam ayrıştırmaları olarak modelleyen doğrusal olmayan ARDL (NARDL) çerçevesini önermektedir.

Doğrusal olmayan ilişkiler, dağıtılmış gecikme değişkenlerinin pozitif ve negatif kısmi toplam ayrıştırmaları olarak modellenir. Herhangi bir verilen bir \(z_t\) değişkeninin \(z_t = z_0 + z_t^+ + z_t^-\) şeklinde ayrıştırılabileceğini hatırlayın.

\[ z_t^+ = \sum_{s=1}^{t} \max(\Delta z_s, 0) \] \[ z_t^- = \sum_{s=1}^{t} \min(\Delta z_s, 0) \] \(z_t\) bir dağıtılmış gecikme değişkeni olduğunda ve \(z_t^+ \neq z_t^-\) durumunda, dağıtılmış gecikme değişkeni asimetrik etkiler gösterir; bu durumda pozitif değişimlerin bağımlı değişken üzerindeki etkisi, negatif karşılıklarından farklıdır. Diğer yandan,\(z_t^+ = z_t^-\) olduğunda, dağıtılmış gecikme değişkeni bağımlı değişken üzerinde simetrik etkiler gösterir ve klasik ARDL etkisine indirgenir.

NARDL çerçevesi ayrıca asimetrik dinamik çarpanlar sunar. VAR literatüründeki etki-tepki eğrilerine benzer olan bu yapılar, her bir doğrusal olmayan dağıtılmış gecikme regresörünün uzun dönemli (eşbütünleşme) durumuna uyumunun asimetrik yollarını izler.

Aşağıda, NARDL çerçevesini Bosna’nın turizm sektörünü ülke ekonomisinin durumuna bağlayan uzun dönemli (eşbütünleşme) ve kısa dönemli (uyarlama) dinamikleri tanımlamak için uygulayacağız.

Veri ve Motivasyon

Analizi gerçekleştirmek için verileri doğrudan Bosna Hersek İstatistik Ajansı’ndan (ASBH) toplayacağız. Özellikle 5 farklı zaman serisiyle ilgileniyoruz:

GDP: gayri safi yurt içi hasıla FTA: yabancı turist gelişleri, foreign tourist arrivals. DTA: yerli turist gelişleri, domestic tourist arrivals FTS: yabancı turist konaklama süresi, foreign tourist length of stay DTS: yerli turist konaklama süresi, domestic tourist length of stay

ASBH’nin yabancı turisti “BH dışında daimi ikameti olan ve BH’de geçici olarak ikamet eden ve bir otelde veya [benzer] konaklama tesisinde en az bir gece geçiren kişi” olarak tanımladığını unutmayın. Benzer şekilde, yerli turisti “BH içinde daimi ikameti olan ve ikamet yeri dışında bir otelde veya [benzer] konaklama tesisinde en az bir gece geçiren kişi” olarak tanımlar. Ayrıca turist gelişlerini “bir konaklama tesisine gelen ve konaklama kaydı yapan kişilerin (turistlerin) sayısı” ve turist konaklama süresini “bir konaklama tesisinde bir kişinin (turistin) kayıtlı gecelik konaklamaları” olarak tanımlar.

Verilerimizle ilgili birkaç uyarıdan bahsetmek de önemlidir. İlk olarak, önceki yıl fiyatları kullanılarak harcama yaklaşımıyla ölçülen Bosna’nın GSYİH’si, üç aylık olarak toplanmaktadır. Ayrıca, GSYİH verileri 2000 1. çeyrekten 2021 4. çeyreğe kadar mevcut olsa da, ASBH 2021’den sonraki verileri tahmin edilen ve gerçek olmayan olarak etiketlemiştir. Bir önlem olarak, bu verileri tamamen göz ardı etmeye ve GSYİH serisini 2020’nin 4. çeyreğinde yapılan son gerçek ölçümleri kapsayacak şekilde kısaltmaya karar verdik.

Buna karşılık, turizm sektörü değişkenleri aylık olarak toplanmaktadır ve Ocak 2008’e kadar uzanmaktadır. Bu, tüm değişkenleri aynı çerçevede eş zamanlı olarak kullanma açısından açıkça bir zorluk teşkil etmektedir. Özellikle, turizm değişkenleri arasında mevcut olan daha uzun zaman serilerinden yararlanmak için, GSYİH’yi düşük üç aylık frekansından Denton yöntemini kullanarak daha yüksek aylık frekansa dönüştürmeye karar verdik. Kolaylık olması açısından, bu verilerin önceden işlenmiş bir versiyonunu buradan indirilebilecek bir EViews çalışma dosyası olarak sunuyoruz.

Herhangi bir ileri analiz yapmadan önce, ele aldığımız verileri anlamak da teşvik edilir. Bu, sadece Bosna’nın turizm sektörünün durumu hakkında bir fikir vermekle kalmayacak, aynı zamanda daha sonra kullanmaya çalışabileceğimiz anlamlı modelleri belirlememize de yardımcı olacaktır.

Bosna’nın yabancı ve yerli turist gelişlerinin her ikisi de 2008 ve 2022 arasındaki raporlama döneminde önemli bir hareketlilik görmüş olsa da, yabancı turizm neredeyse üstel bir dönüşüm geçirmiş görünmektedir. Bu belirgin kontrast aşağıdaki Şekil 1a’da gösterilmiştir. Aslında (Şekil 1b’ye bakın), COVID-19 pandemi yılları öncesindeki dönemde, yıllık turist gelişlerinin ortalama yıllık yüzde değişimi sırasıyla yabancı ve yerli sektörler için %13 ve %4,12 idi. 2015 yılında, turist gelişlerindeki yıllık büyüme yabancılar arasında %26,53, yerliler arasında %13,1 civarındaydı.

bhgdp <- read_csv("bhgdp.csv")
bh <- read_csv("bhdata.csv")

gdp <- bhgdp[33:84,]

library(tempdisagg)

## Warning: package 'tempdisagg' was built under R version 4.4.3

library(dplyr)
library(zoo)

# Convert your data to time series format
# First, convert Date to proper format
gdp$Date <- as.yearqtr(gdp$Date, format = "%m/%d/%Y")

# Create a ts object with frequency = 4 (quarterly)
gdp_ts <- ts(gdp$gdp, 
             start = 2008, 
             frequency = 4)

# Use the Denton method for temporal disaggregation
# Without an indicator, we'll use the proportional Denton method
gdp_monthly <- td(gdp_ts ~ 1, to = 12, method = "denton-cholette", criterion = "proportional")

# Extract the results
monthly_gdp <- predict(gdp_monthly)

# Convert back to a data frame with dates
start_date <- as.Date(as.yearmon(gdp$Date[1]))
end_date <- as.Date(as.yearmon(gdp$Date[nrow(gdp)]) + 2/12) # End of the last quarter

monthly_dates <- seq.Date(from = start_date, to = end_date, by = "month")

monthly_gdp_df <- data.frame(
  Date = monthly_dates,
  gdp = as.numeric(monthly_gdp)
)

bdf <- bh[1:156,]

bdf$gdp <- monthly_gdp_df$gdp

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

bdf$Year <- substr(bdf$EN, 1, 4)

# Aggregate TA and DA by year (sum values for each year)
yearly_data <- bdf %>%
  group_by(Year) %>%
  summarize(
    TA_Total = sum(TA, na.rm = TRUE),
    DA_Total = sum(DA, na.rm = TRUE)
  )

# Convert data to long format for easier plotting
library(tidyr)
yearly_data_long <- yearly_data %>%
  pivot_longer(
    cols = c(TA_Total, DA_Total),
    names_to = "Variable",
    values_to = "Value"
  )

# Create a bar graph
ggplot(yearly_data_long, aes(x = Year, y = Value, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(
    values = c("TA_Total" = "#3366CC", "DA_Total" = "#FF9933"),
    labels = c("TA_Total" = "TA", "DA_Total" = "DA")
  ) +
  labs(
    title = "Yearly Comparison of TA and DA Values",
    x = "Year",
    y = "Total Value",
    fill = "Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "top"
  ) +
  scale_y_continuous(labels = scales::comma)

yearly_pct_change <- yearly_data %>%
  arrange(Year) %>%
  mutate(
    TA_Pct_Change = (TA_Total / lag(TA_Total) - 1) * 100,
    DA_Pct_Change = (DA_Total / lag(DA_Total) - 1) * 100
  ) %>%
  # Remove the first row which will have NA values
  filter(!is.na(TA_Pct_Change))

# Convert to long format for plotting
yearly_pct_change_long <- yearly_pct_change %>%
  select(Year, TA_Pct_Change, DA_Pct_Change) %>%
  pivot_longer(
    cols = c(TA_Pct_Change, DA_Pct_Change),
    names_to = "Variable",
    values_to = "Percent_Change"
  )

# Modify the variable names for better display
yearly_pct_change_long$Variable <- gsub("_Pct_Change", "", yearly_pct_change_long$Variable)

# Create the bar graph
ggplot(yearly_pct_change_long, aes(x = Year, y = Percent_Change, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_fill_manual(
    values = c("TA" = "#3366CC", "DA" = "#FF9933")
  ) +
  labs(
    title = "Year-on-Year Percentage Change in TA and DA",
    subtitle = "Comparison of annual growth rates",
    x = "Year",
    y = "Percentage Change (%)",
    fill = "Variable"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "top"
  ) +
  scale_y_continuous(labels = function(x) paste0(round(x, 1), "%"))

# Convert Date column to proper date format
bdf$Dt <- as.Date(paste0("01-", bdf$Date), format = "%d-%b-%y")

# Create a time series plot
ggplot(bdf, aes(x = Dt)) +
  geom_line(aes(y = TA, color = "TA"), linewidth = 1) +
  geom_line(aes(y = DA, color = "DA"), linewidth = 1) +
  scale_color_manual(
    values = c("TA" = "#3366CC", "DA" = "#FF9933"),
    name = "Variable"
  ) +
  labs(
    title = "Time Series of TA and DA Values",
    subtitle = "Monthly data from 2008 onwards",
    x = "Date",
    y = "Value"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
  scale_y_continuous(labels = scales::comma)

Şekillerde göze çarpan şey, hem yerli hem de yabancı turist varışlarının öngörülebilir mevsimsel etkiler göstermesidir. Her ikisi de yılın başlarında dip noktalara ulaşır ve sırasıyla turizm sezonunun başında ve sonunda zirvelere ulaşır. Daha da önemlisi, yabancı turistler arasında zirve turizm aylarındaki varışların büyüklüğü, yerli eşdeğerinin çok üzerindedir.

bdf <- bdf %>% mutate(mtn = TN/TA,
                      mfn = FN/FA,
                      mdn = DN/DA)

library(explore)
bdf %>% describe_all()

## # A tibble: 14 × 8
##    variable type     na na_pct unique        min       mean        max
##    <chr>    <chr> <int>  <dbl>  <int>      <dbl>      <dbl>      <dbl>
##  1 EN       chr       0      0    156      NA         NA         NA   
##  2 TA       dbl       0      0    156    1539      77302.    204601   
##  3 DA       dbl       0      0    155    1271      27802.     52762   
##  4 FA       dbl       0      0    156     268      49500.    165554   
##  5 Date     chr       0      0    156      NA         NA         NA   
##  6 TN       dbl       0      0    156    7851     164053.    454953   
##  7 DN       dbl       0      0    156    5541      60265.    158048   
##  8 FN       dbl       0      0    156    1670     103795.    357593   
##  9 gdp      dbl       0      0    156 1821123.   2495816.   3117643.  
## 10 Year     chr       0      0     13      NA         NA         NA   
## 11 Dt       dat       0      0    156      NA         NA         NA   
## 12 mtn      dbl       0      0    156       1.79       2.16       5.1 
## 13 mfn      dbl       0      0    156       1.77       2.2        8.62
## 14 mdn      dbl       0      0    156       1.66       2.18       4.36

mms <- bdf %>% filter(Dt<"2020-01-01") %>%
  summarize(
    mdn_mean = mean(mdn, na.rm = TRUE),
    mdn_median = median(mdn, na.rm = TRUE),
    mdn_sd = sd(mdn, na.rm = TRUE),
    mfn_mean = mean(mfn, na.rm = TRUE),
    mfn_median = median(mfn, na.rm = TRUE),
    mfn_sd = sd(mfn, na.rm = TRUE)
  )

mms

## # A tibble: 1 × 6
##   mdn_mean mdn_median mdn_sd mfn_mean mfn_median mfn_sd
##      <dbl>      <dbl>  <dbl>    <dbl>      <dbl>  <dbl>
## 1     2.15       2.09  0.328     2.13       2.09  0.210

# First, let's extract the month from the Dt column
bdf$month <- month(bdf$Dt)

# Filter for Q1 months (January, February, March) and years before 2020
q1_before_2020 <- bdf %>%
  filter(month %in% c(1, 2, 3) & Year < "2020")

q3_before_2020 <- bdf %>%
  filter(month %in% c(7, 8, 9) & Year < "2020")

# Calculate the mean and median for mdn and mfn
result <- q1_before_2020 %>%
  summarize(
    mdn_mean = mean(mdn, na.rm = TRUE),
    mdn_median = median(mdn, na.rm = TRUE),
    mfn_mean = mean(mfn, na.rm = TRUE),
    mfn_median = median(mfn, na.rm = TRUE)
  )

result2 <- q3_before_2020 %>%
  summarize(
    mdn_mean = mean(mdn, na.rm = TRUE),
    mdn_median = median(mdn, na.rm = TRUE),
    mfn_mean = mean(mfn, na.rm = TRUE),
    mfn_median = median(mfn, na.rm = TRUE)
  )

result

## # A tibble: 1 × 4
##   mdn_mean mdn_median mfn_mean mfn_median
##      <dbl>      <dbl>    <dbl>      <dbl>
## 1     1.98       1.95     2.33       2.35

result2

## # A tibble: 1 × 4
##   mdn_mean mdn_median mfn_mean mfn_median
##      <dbl>      <dbl>    <dbl>      <dbl>
## 1     2.56       2.52     2.20       2.19

Turizm sıklıkla takdirî zamana sıkı sıkıya bağlı olduğundan, çalışan nüfusun büyük bir kısmı için bu kaynak uzun yıllar boyunca sabit kalacaktır. Bu nedenle, Bosna’ya yapılan ziyaretlerin uzunluğunu belirleyen davranış kalıplarının raporlama dönemimiz boyunca çok az değişmesini bekliyoruz. Özellikle, COVID-19 salgınından önceki dönemde, ortalama (medyan) kalış süresi yerli ve yabancı turistler için sırasıyla 2,15 (2,08) ve 2,13 (2,08) gündü, standart sapmalar ise sırasıyla 0,32 ve 0,21 civarında seyretti. Şekil 5b’de, yerli ve yabancı turist konaklamalarının mevsimselliğini anlamaya çalışıyoruz. Turist varışlarında olduğu gibi, bu eğriler de mevsimsellik göstermektedir. Bununla birlikte, yerli turistlerin kış aylarına kıyasla yaz aylarında daha uzun süre kalırken, yabancı turistler için bunun tam tersinin geçerli olduğunu gösteren ilginç bir kalıp bulunmaktadır. Aslında 2020 öncesi yıllarda, 1. çeyrekte yerli ve yabancı turistler için ortalama (medyan) kalış süresi sırasıyla 1,98 (1,95) ve 2,33 (2,35) gündü. Öte yandan, aynı dönemde, 3. çeyrekte yerli ve yabancı turistler için ortalama (medyan) kalış süresi sırasıyla 2,56 (2,52) ve 2,20 (2,19) gündü.

Şimdi Bosna’nın GSYİH’sine bir göz atalım. Aşağıdaki Şekil, yabancı turist varışlarının etkilerini yakından taklit eden mevsimsel etkilerle trend büyümeyi belirleyebiliriz; yani, GSYİH Ocak ayında azalır ve Temmuz ayında zirve yapar. Daha fazla bilgi için, diğer Şekil standartlaştırılmış aylık GSYİH’yi ve hem yerli hem de yabancı turist varışlarını çizer. Bu çizimin gösterdiği şey, turist varışlarının ortalama değerleri etrafında GSYİH’den önemli ölçüde daha fazla dalgalanmasıdır. Dahası, yerel turizm raporlanan örneğin başında bu dalgalanmalara hakim gibi görünürken, yabancı turizmdeki sapmalar raporlanan örneğin sonuna hakim gibi görünmektedir. Daha da önemlisi, GSYİH ile yabancı turist varışları arasında güçlü bir pozitif korelasyon görüyoruz.

ggplot(bdf, aes(Dt, gdp)) + geom_line()

library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:readr':
## 
##     col_factor

# Create a function to standardize variables
# This converts variables to z-scores: (value - mean) / sd
standardize <- function(x) {
  return((x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE))
}

# Create a new dataframe with standardized variables
std_df <- bdf %>%
  mutate(
    gdp_std = standardize(gdp),
    mdn_std = standardize(DA),
    mfn_std = standardize(FA)
  )

# Convert to long format for easier plotting
std_long <- std_df %>%
  select(Dt, gdp_std, mdn_std, mfn_std) %>%
  pivot_longer(
    cols = c(gdp_std, mdn_std, mfn_std),
    names_to = "Variable",
    values_to = "Standardized_Value"
  )

# Clean up variable names for the legend
std_long$Variable <- case_when(
  std_long$Variable == "gdp_std" ~ "GDP",
  std_long$Variable == "mdn_std" ~ "MDN",
  std_long$Variable == "mfn_std" ~ "MFN",
  TRUE ~ std_long$Variable
)

# Create the plot
ggplot(std_long, aes(x = Dt, y = Standardized_Value, color = Variable, group = Variable)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_color_manual(
    values = c("GDP" = "#1B9E77", "MDN" = "#D95F02", "MFN" = "#7570B3")
  ) +
  labs(
    title = "Standardized Comparison of GDP, MDN, and MFN Over Time",
    subtitle = "Variables converted to z-scores for comparison",
    x = "Date",
    y = "Standardized Value (z-score)",
    color = "Variable"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(labels = function(x) sprintf("%.1f", x))

# Optional: Create a correlation matrix to show relationships
cor_matrix <- cor(std_df %>% select(gdp_std, mdn_std, mfn_std), use = "complete.obs")
print(cor_matrix)

##           gdp_std   mdn_std   mfn_std
## gdp_std 1.0000000 0.5851591 0.5756952
## mdn_std 0.5851591 1.0000000 0.6476426
## mfn_std 0.5756952 0.6476426 1.0000000

Tahmine geçmeden önce, serilerimizin entegrasyon sıralarını incelemek de ihtiyatlı olacaktır. Girişte belirtildiği gibi, ARDL tahmini I(2) değişkenlerinin varlığında geçerli değildir, ancak I(0) ve I(1) değişkenlerinin bir karışımını barındırır. İlk test, seviyelerdeki serilerde sabit ve trend ile Im, Pesaran ve Shin (IPS) testini gerçekleştirecektir. Bu testin sıfır hipotezi bir birim kök olduğundan, sıfıra yakın p-değerlerine sahip testler sıfırı reddedecek ve seriyi I(0) olarak tanımlayacaktır. Aksi takdirde, sıfırdan önemli ölçüde uzakta p-değerlerine sahip testler sıfırı reddetmekte ve seriyi I(1) olarak tanımlamakta başarısız olacaktır. Genişletme yoluyla, I(2) serileri serinin ilk farkı üzerindeki testleri tekrarlayarak belirlenebilir.

# Load required packages
library(urca)
library(knitr)

# Function to perform multiple unit root tests on a single series
test_unit_root <- function(series, series_name) {
  # Skip if there are NA values or insufficient data
  if (sum(is.na(series)) > 0 || length(series) < 10) {
    cat("\nSkipping", series_name, "due to NA values or insufficient data\n")
    return(NULL)
  }
  
  # 1. ADF test with constant and trend
  adf_result <- tryCatch({
    ur.df(series, type = "trend", lags = 4, selectlags = "AIC")
  }, error = function(e) {
    cat("\nError in ADF test for", series_name, ":", e$message, "\n")
    return(NULL)
  })
  
  if (is.null(adf_result)) {
    return(NULL)
  }
  
  # Extract test statistic and critical values
  adf_test_stat <- adf_result@teststat[1]
  adf_crit_vals <- adf_result@cval[1,]
  
  # 2. Phillips-Perron test with constant and trend
  pp_result <- tryCatch({
    ur.pp(series, type = "Z-tau", model = "trend", lags = "short")
  }, error = function(e) {
    cat("\nError in PP test for", series_name, ":", e$message, "\n")
    return(NULL)
  })
  
  if (is.null(pp_result)) {
    pp_test_stat <- NA
    pp_crit_vals <- c(NA, NA, NA)
  } else {
    pp_test_stat <- pp_result@teststat
    pp_crit_vals <- pp_result@cval
  }
  
  # 3. KPSS test with trend (null: stationary)
  kpss_result <- tryCatch({
    kpss.test(series, null = "Trend")
  }, error = function(e) {
    cat("\nError in KPSS test for", series_name, ":", e$message, "\n")
    return(NULL)
  })
  
  if (is.null(kpss_result)) {
    kpss_stat <- NA
    kpss_p_val <- NA
  } else {
    kpss_stat <- kpss_result$statistic
    kpss_p_val <- kpss_result$p.value
  }
  
  # Return the test results in a list
  return(list(
    series_name = series_name,
    adf = list(statistic = adf_test_stat, critical_values = adf_crit_vals),
    pp = list(statistic = pp_test_stat, critical_values = pp_crit_vals),
    kpss = list(statistic = kpss_stat, p_value = kpss_p_val)
  ))
}

# Create a function to process both level and differenced series
process_series <- function(data, series_names) {
  # Create results dataframes
  level_results <- data.frame(
    Series = character(),
    ADF_Statistic = numeric(),
    ADF_5pct_Critical = numeric(),
    ADF_Result = character(),
    PP_Statistic = numeric(),
    PP_5pct_Critical = numeric(),
    PP_Result = character(),
    KPSS_Statistic = numeric(),
    KPSS_p_value = numeric(),
    KPSS_Result = character(),
    Overall = character(),
    stringsAsFactors = FALSE
  )
  
  diff_results <- level_results
  
  # Process each series
  for (series_name in series_names) {
    cat("\n===== Processing", series_name, "=====\n")
    
    # Get the level series data
    level_data <- data[[series_name]]
    
    # Run tests on level data
    cat("\nTesting level data for", series_name, "\n")
    level_test_results <- test_unit_root(level_data, series_name)
    
    if (!is.null(level_test_results)) {
      # Determine results
      adf_result <- ifelse(level_test_results$adf$statistic < level_test_results$adf$critical_values[2], 
                          "Stationary", "Non-stationary")
      
      pp_result <- ifelse(level_test_results$pp$statistic < level_test_results$pp$critical_values[2], 
                         "Stationary", "Non-stationary")
      
      kpss_result <- ifelse(level_test_results$kpss$p_value > 0.05, 
                           "Stationary", "Non-stationary")
      
      # Determine overall result
      if ((adf_result == "Non-stationary" || pp_result == "Non-stationary") && 
          kpss_result == "Non-stationary") {
        overall <- "Non-stationary"
      } else if (adf_result == "Stationary" && pp_result == "Stationary" && 
                kpss_result == "Stationary") {
        overall <- "Stationary"
      } else {
        overall <- "Mixed results"
      }
      
      # Add to level results dataframe
      level_results <- rbind(level_results, data.frame(
        Series = series_name,
        ADF_Statistic = level_test_results$adf$statistic,
        ADF_5pct_Critical = level_test_results$adf$critical_values[2],
        ADF_Result = adf_result,
        PP_Statistic = level_test_results$pp$statistic,
        PP_5pct_Critical = level_test_results$pp$critical_values[2],
        PP_Result = pp_result,
        KPSS_Statistic = level_test_results$kpss$statistic,
        KPSS_p_value = level_test_results$kpss$p_value,
        KPSS_Result = kpss_result,
        Overall = overall,
        stringsAsFactors = FALSE
      ))
    }
    
    # Create differenced series
    diff_data <- diff(level_data, lag = 1, differences = 1)
    diff_data <- diff_data[!is.na(diff_data)]
    
    # Run tests on differenced data
    cat("\nTesting first difference for", series_name, "\n")
    diff_test_results <- test_unit_root(diff_data, paste0("d_", series_name))
    
    if (!is.null(diff_test_results)) {
      # Determine results
      adf_result <- ifelse(diff_test_results$adf$statistic < diff_test_results$adf$critical_values[2], 
                          "Stationary", "Non-stationary")
      
      pp_result <- ifelse(diff_test_results$pp$statistic < diff_test_results$pp$critical_values[2], 
                         "Stationary", "Non-stationary")
      
      kpss_result <- ifelse(diff_test_results$kpss$p_value > 0.05, 
                           "Stationary", "Non-stationary")
      
      # Determine overall result
      if ((adf_result == "Stationary" || pp_result == "Stationary") && 
          kpss_result == "Stationary") {
        overall <- "Stationary"
      } else if (adf_result == "Non-stationary" && pp_result == "Non-stationary" && 
                kpss_result == "Non-stationary") {
        overall <- "Non-stationary"
      } else {
        overall <- "Mixed results"
      }
      
      # Add to differenced results dataframe
      diff_results <- rbind(diff_results, data.frame(
        Series = paste0("d_", series_name),
        ADF_Statistic = diff_test_results$adf$statistic,
        ADF_5pct_Critical = diff_test_results$adf$critical_values[2],
        ADF_Result = adf_result,
        PP_Statistic = diff_test_results$pp$statistic,
        PP_5pct_Critical = diff_test_results$pp$critical_values[2],
        PP_Result = pp_result,
        KPSS_Statistic = diff_test_results$kpss$statistic,
        KPSS_p_value = diff_test_results$kpss$p_value,
        KPSS_Result = kpss_result,
        Overall = overall,
        stringsAsFactors = FALSE
      ))
    }
  }
  
  return(list(level_results = level_results, diff_results = diff_results))
}

# Variables to test
series_to_test <- c("DA", "FA", "DN", "FN", "gdp", "mfn", "mdn")

# Process all series
results <- process_series(bdf, series_to_test)

## 
## ===== Processing DA =====
## 
## Testing level data for DA

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## Testing first difference for DA

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing FA =====
## 
## Testing level data for FA

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## Testing first difference for FA

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing DN =====
## 
## Testing level data for DN

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## Testing first difference for DN

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing FN =====
## 
## Testing level data for FN

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## Testing first difference for FN

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing gdp =====
## 
## Testing level data for gdp

## Warning in kpss.test(series, null = "Trend"): p-value smaller than printed
## p-value

## 
## Testing first difference for gdp

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing mfn =====
## 
## Testing level data for mfn 
## 
## Testing first difference for mfn

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

## 
## ===== Processing mdn =====
## 
## Testing level data for mdn 
## 
## Testing first difference for mdn

## Warning in kpss.test(series, null = "Trend"): p-value greater than printed
## p-value

# Print summary tables
cat("\n===== Summary of Unit Root Tests for Level Series =====\n")

## 
## ===== Summary of Unit Root Tests for Level Series =====

print(kable(results$level_results))

## 
## 
## |      |Series | ADF_Statistic| ADF_5pct_Critical|ADF_Result | PP_Statistic| PP_5pct_Critical|PP_Result  | KPSS_Statistic| KPSS_p_value|KPSS_Result    |Overall       |
## |:-----|:------|-------------:|-----------------:|:----------|------------:|----------------:|:----------|--------------:|------------:|:--------------|:-------------|
## |5pct  |DA     |     -4.979994|             -3.43|Stationary |    -6.466957|          -3.4394|Stationary |      0.0690760|    0.1000000|Stationary     |Stationary    |
## |5pct1 |FA     |     -5.938371|             -3.43|Stationary |    -3.999862|          -3.4394|Stationary |      0.1055322|    0.1000000|Stationary     |Stationary    |
## |5pct2 |DN     |     -8.357627|             -3.43|Stationary |    -5.518290|          -3.4394|Stationary |      0.0467794|    0.1000000|Stationary     |Stationary    |
## |5pct3 |FN     |     -5.572218|             -3.43|Stationary |    -4.151980|          -3.4394|Stationary |      0.1098142|    0.1000000|Stationary     |Stationary    |
## |5pct4 |gdp    |     -5.313247|             -3.43|Stationary |    -4.248710|          -3.4394|Stationary |      0.2936156|    0.0100000|Non-stationary |Mixed results |
## |5pct5 |mfn    |     -7.825397|             -3.43|Stationary |    -9.333400|          -3.4394|Stationary |      0.1380141|    0.0647888|Stationary     |Stationary    |
## |5pct6 |mdn    |     -7.027968|             -3.43|Stationary |    -6.074099|          -3.4394|Stationary |      0.1743774|    0.0263522|Non-stationary |Mixed results |

cat("\n===== Summary of Unit Root Tests for First Differences =====\n")

## 
## ===== Summary of Unit Root Tests for First Differences =====

print(kable(results$diff_results))

## 
## 
## |      |Series | ADF_Statistic| ADF_5pct_Critical|ADF_Result | PP_Statistic| PP_5pct_Critical|PP_Result  | KPSS_Statistic| KPSS_p_value|KPSS_Result |Overall    |
## |:-----|:------|-------------:|-----------------:|:----------|------------:|----------------:|:----------|--------------:|------------:|:-----------|:----------|
## |5pct  |d_DA   |     -9.626264|             -3.43|Stationary |   -14.067250|        -3.439579|Stationary |      0.0159488|          0.1|Stationary  |Stationary |
## |5pct1 |d_FA   |     -7.926351|             -3.43|Stationary |    -9.193746|        -3.439579|Stationary |      0.0210457|          0.1|Stationary  |Stationary |
## |5pct2 |d_DN   |     -9.021012|             -3.43|Stationary |   -10.143162|        -3.439579|Stationary |      0.0140405|          0.1|Stationary  |Stationary |
## |5pct3 |d_FN   |     -7.494099|             -3.43|Stationary |   -10.231043|        -3.439579|Stationary |      0.0218592|          0.1|Stationary  |Stationary |
## |5pct4 |d_gdp  |     -9.517332|             -3.43|Stationary |    -5.187582|        -3.439579|Stationary |      0.0340251|          0.1|Stationary  |Stationary |
## |5pct5 |d_mfn  |     -9.031583|             -3.43|Stationary |   -21.581164|        -3.439579|Stationary |      0.0179140|          0.1|Stationary  |Stationary |
## |5pct6 |d_mdn  |     -8.793459|             -3.43|Stationary |   -12.757103|        -3.439579|Stationary |      0.0152485|          0.1|Stationary  |Stationary |

# Determine the order of integration
integration_order <- data.frame(
  Series = series_to_test,
  Level_Stationary = results$level_results$Overall == "Stationary",
  Diff_Stationary = NA,
  Order_of_Integration = NA,
  stringsAsFactors = FALSE
)

for (i in 1:nrow(integration_order)) {
  series_name <- integration_order$Series[i]
  diff_row <- which(results$diff_results$Series == paste0("d_", series_name))
  
  if (length(diff_row) > 0) {
    integration_order$Diff_Stationary[i] <- results$diff_results$Overall[diff_row] == "Stationary"
    
    if (integration_order$Level_Stationary[i]) {
      integration_order$Order_of_Integration[i] <- "I(0)"
    } else if (integration_order$Diff_Stationary[i]) {
      integration_order$Order_of_Integration[i] <- "I(1)"
    } else {
      integration_order$Order_of_Integration[i] <- "Higher than I(1)"
    }
  }
}

cat("\n===== Order of Integration =====\n")

## 
## ===== Order of Integration =====

print(kable(integration_order))

## 
## 
## |Series |Level_Stationary |Diff_Stationary |Order_of_Integration |
## |:------|:----------------|:---------------|:--------------------|
## |DA     |TRUE             |TRUE            |I(0)                 |
## |FA     |TRUE             |TRUE            |I(0)                 |
## |DN     |TRUE             |TRUE            |I(0)                 |
## |FN     |TRUE             |TRUE            |I(0)                 |
## |gdp    |FALSE            |TRUE            |I(1)                 |
## |mfn    |TRUE             |TRUE            |I(0)                 |
## |mdn    |FALSE            |TRUE            |I(1)                 |

Seviyelerdeki birim kök testleri, yerli turistlerin kalış sürelerini ve GSYİH’yi 1. dereceden entegre olarak tanımlarken, kalan değişkenler 0. dereceden entegre olarak tanımlanmaktadır. Öte yandan, birinci farklardaki birim kök testleri için her sıfır hipotezi reddedilir ve bu nedenle hiçbir değişken I(2) olarak tanımlanmaz.

Asimetri Tahmini

Bu bölümdeki amacımız, Bosna’nın turist değişkenleri DTA, DTS, FTA, FTS’nin, değişken GSYİH olarak ölçülen ekonominin durumuyla ilişkisini bağlayan bir NARDL modeli tahmin etmektir. Resmi olarak, aşağıda koşullu hata düzeltme (CEC) biçiminde ifade edilen NARDL modeli şudur:

\(\Delta \ln(GDP) = \phi_{\text{GDP}} \ln(\text{GDP})_{t-1}\) \(+\phi_{\text{DTA}}^{+} \ln(\text{DTA})_{t-1}^{+} + \phi_{\text{DTA}}^{-} \ln(\text{DTA})_{t-1}^{-} + \phi_{\text{FTA}}^{+} \ln(\text{FTA})_{t-1}^{+} + \phi_{\text{FTA}}^{-} \ln(\text{FTA})_{t-1}^{-}\) \(+\phi_{\text{DTS}}^{+}\text{DTS}_{t-1}^{+} + \phi_{\text{DTS}}^{-}\text{DTS}_{t-1}^{-} + \phi_{\text{FTS}}^{+}\text{FTS}_{t-1}^{+} + \phi_{\text{FTS}}^{-}\text{FTS}_{t-1}^{-}\) \(+ \sum_{j=1}^{p-1} \gamma_{\text{GDP}, j}\Delta \ln(\text{GDP})_{t-j}\) \(+ \sum_{k_1=1}^{q_1-1} \left( \gamma_{\text{DTA}, k_1}^{+} \Delta \ln(\text{DTA})_{t-k_1}^{+} + \gamma_{\text{DTA}, k_1}^{-} \Delta \ln(\text{DTA})_{t-k_1}^{-} \right) + \sum_{k_2=1}^{q_2-1} \left( \gamma_{\text{FTA}, k_2}^{+} \Delta \ln(\text{FTA})_{t-k_2}^{+} + \gamma_{\text{FTA}, k_2}^{-} \Delta \ln(\text{FTA})_{t-k_2}^{-} \right)\) \(+ \sum_{k_3=1}^{q_3-1} \left( \gamma_{\text{DTS}, k_3}^{+} \Delta \ln(\text{DTS})_{t-k_3}^{+} + \gamma_{\text{DTS}, k_3}^{-} \Delta \ln(\text{DTS})_{t-k_3}^{-} \right) + \sum_{k_4=1}^{q_4-1} \left( \gamma_{\text{FTS}, k_4}^{+} \Delta \ln(\text{FTS})_{t-k_4}^{+} + \gamma_{\text{FTS}, k_4}^{-} \Delta \ln(\text{FTS})_{t-k_4}^{-} \right)\) \(+\alpha_0 + \alpha_1 t + \sum_{i=1}^{11} \delta_i m_i + \epsilon_t\)

Bu, GSYİH’nin p mertebesinde bir otoregresif süreç olarak girdiği ve DTA, FTA, DTS, FTS’nin sırasıyla q1, q2, q3, q4 mertebelerinde asimetrik dağıtılmış gecikmeli değişkenler olarak girdiği bir NARDL modelini açıklar. İkincisi, sırasıyla eşbütünleşmeyi (seviyeler veya uzun vadeli), ayarlamayı (farklar veya kısa vadeli) ve deterministik (mevsimsellik) dinamiklerini karakterize eder. Ayrıca, + ve - üst simgelere sahip değişkenler, sırasıyla, altta yatan dağıtılmış gecikmeli değişkenin pozitif ve negatif kısmi toplam ayrıştırmalarını belirtir. Buradaki pozitif ve negatif kısmi toplamlar, Bosna’nın turizm sektöründeki asimetrilerin hem uzun hem de kısa vadede ekonomik gelişimine nasıl yansıdığını açıkça modellemektedir.

Deterministik dinamikler de kısa bir yorumu hak ediyor. Özellikle, α0 ve α1 sırasıyla sabit ve doğrusal eğilimin etkisini yakalar. Aylık mevsimselliği de yakalamak için, katsayılar δi, ay i için mevsimsel kukla değişkenler mi ile ilişkilidir.

Analize yukarıdaki modeli tahmin ederek başlayacağız ve tüm turizm değişkenlerini hem ayarlayan hem de eşbütünleşik dinamiklerde asimetrik olarak ele alacağız. Otoregresif ve dağıtılmış gecikme sıralarını belirlemek için, bağımlı değişken ve regresörlerin her biri için en fazla 3 gecikmeye izin veren otomatik gecikme seçimi gerçekleştireceğiz (varsayılan seçenekler). Bu, tüm değişkenlerin geçmişte en fazla 3 dönemdeki (ay) değerlere; başka bir deyişle, tek bir çeyreğe bağlı olduğunu etkili bir şekilde belirtir.

Pandemi yıllarının karmaşıklıklarından kaçınmak için, yalnızca 2020’den önceki yıllara ait verileri kullanacağız.

# Gerekli kütüphaneleri yükleyelim
library(dplyr)
library(nardl) # NARDL analizi için

## Warning: package 'nardl' was built under R version 4.4.3

# Logaritmik dönüşümleri yapalım
bdf$log_DA <- log(bdf$DA)
bdf$log_FA <- log(bdf$FA)
bdf$log_gdp <- log(bdf$gdp)

# Değişkenlerin pozitif ve negatif bileşenlerini oluşturalım
# NARDL için parçalama fonksiyonu
decompose_var <- function(x) {
  # Farkları hesapla
  dx <- c(0, diff(x))
  
  # Pozitif ve negatif değişimleri ayır
  pos_changes <- pmax(dx, 0)
  neg_changes <- pmin(dx, 0)
  
  # Kümülatif toplamları hesapla
  pos_cum <- cumsum(pos_changes)
  neg_cum <- cumsum(neg_changes)
  
  return(list(positive = pos_cum, negative = neg_cum))
}

# log(DA)'nın pozitif ve negatif bileşenleri
da_decomp <- decompose_var(bdf$log_DA)
bdf$log_DA_pos <- da_decomp$positive
bdf$log_DA_neg <- da_decomp$negative

# log(FA)'nın pozitif ve negatif bileşenleri
fa_decomp <- decompose_var(bdf$log_FA)
bdf$log_FA_pos <- fa_decomp$positive
bdf$log_FA_neg <- fa_decomp$negative

# mdn'in pozitif ve negatif bileşenleri
mdn_decomp <- decompose_var(bdf$mdn)
bdf$mdn_pos <- mdn_decomp$positive
bdf$mdn_neg <- mdn_decomp$negative

# mfn'in pozitif ve negatif bileşenleri
mfn_decomp <- decompose_var(bdf$mfn)
bdf$mfn_pos <- mfn_decomp$positive
bdf$mfn_neg <- mfn_decomp$negative

# Sonuçları kontrol edelim
head(bdf[, c("log_DA", "log_DA_pos", "log_DA_neg", 
              "log_FA", "log_FA_pos", "log_FA_neg",
              "mdn", "mdn_pos", "mdn_neg",
              "mfn", "mfn_pos", "mfn_neg")])

## # A tibble: 6 × 12
##   log_DA log_DA_pos log_DA_neg log_FA log_FA_pos log_FA_neg   mdn mdn_pos
##    <dbl>      <dbl>      <dbl>  <dbl>      <dbl>      <dbl> <dbl>   <dbl>
## 1   9.79     0          0        9.81      0         0       2.21   0    
## 2   9.86     0.0698     0        9.79      0        -0.0230  2.01   0    
## 3  10.0      0.220      0        9.91      0.118    -0.0230  1.86   0    
## 4  10.0      0.255      0       10.2       0.395    -0.0230  2.05   0.192
## 5  10.4      0.655      0       10.6       0.765    -0.0230  2.25   0.388
## 6  10.4      0.655     -0.0334  10.4       0.765    -0.172   2.44   0.577
## # ℹ 4 more variables: mdn_neg <dbl>, mfn <dbl>, mfn_pos <dbl>, mfn_neg <dbl>

# 11 kukla değişken oluştur (12. ay referans kategori)
for (i in 1:11) {
  bdf[paste0("month_", i)] <- ifelse(bdf$month == i, 1, 0)
}

bdf$trend <- 1:nrow(bdf)

modnardl <- auto_ardl(log_gdp ~ log_FA_pos + log_FA_neg + log_DA_pos + log_DA_neg + mdn_pos + mdn_neg + mfn_pos + mfn_neg + trend  | month_1 + month_2 + month_3 + month_4 + month_5 + month_6 + month_7 + month_8 + month_9 + month_10 + month_11, data=bdf, max_order = 4)

modnardl$best_order

##    log_gdp log_FA_pos log_FA_neg log_DA_pos log_DA_neg    mdn_pos    mdn_neg 
##          3          3          2          2          0          2          3 
##    mfn_pos    mfn_neg      trend 
##          0          3          3

summary(uecm(modnardl$best_model))

## 
## Time series regression with "ts" data:
## Start = 4, End = 156
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0195675 -0.0027164 -0.0001845  0.0026585  0.0160229 
## 
## Coefficients: (3 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.7160688  0.5305743   5.119 1.26e-06 ***
## L(log_gdp, 1)       -0.1869455  0.0363574  -5.142 1.14e-06 ***
## L(log_FA_pos, 1)     0.0291467  0.0066323   4.395 2.50e-05 ***
## L(log_FA_neg, 1)     0.0202263  0.0060034   3.369  0.00103 ** 
## L(log_DA_pos, 1)    -0.0085789  0.0095109  -0.902  0.36895    
## log_DA_neg          -0.0022581  0.0086124  -0.262  0.79364    
## L(mdn_pos, 1)        0.0142223  0.0087555   1.624  0.10706    
## L(mdn_neg, 1)        0.0076689  0.0064967   1.180  0.24029    
## mfn_pos             -0.0009868  0.0043165  -0.229  0.81958    
## L(mfn_neg, 1)       -0.0039047  0.0042602  -0.917  0.36132    
## L(trend, 1)         -0.0011192  0.0006285  -1.781  0.07759 .  
## d(L(log_gdp, 1))     1.0081445  0.0670649  15.032  < 2e-16 ***
## d(L(log_gdp, 2))    -0.4310539  0.0864805  -4.984 2.24e-06 ***
## d(log_FA_pos)        0.0122738  0.0083583   1.468  0.14473    
## d(L(log_FA_pos, 1)) -0.0080686  0.0077821  -1.037  0.30202    
## d(L(log_FA_pos, 2)) -0.0137138  0.0040623  -3.376  0.00101 ** 
## d(log_FA_neg)        0.0102295  0.0081839   1.250  0.21387    
## d(L(log_FA_neg, 1)) -0.0235599  0.0082976  -2.839  0.00535 ** 
## d(log_DA_pos)       -0.0054329  0.0094946  -0.572  0.56830    
## d(L(log_DA_pos, 1))  0.0211490  0.0090929   2.326  0.02179 *  
## d(mdn_pos)           0.0014559  0.0084396   0.173  0.86335    
## d(L(mdn_pos, 1))     0.0132406  0.0090238   1.467  0.14505    
## d(mdn_neg)           0.0168371  0.0098174   1.715  0.08906 .  
## d(L(mdn_neg, 1))     0.0018430  0.0085811   0.215  0.83033    
## d(L(mdn_neg, 2))     0.0129556  0.0079489   1.630  0.10589    
## d(mfn_neg)          -0.0009353  0.0060069  -0.156  0.87654    
## d(L(mfn_neg, 1))    -0.0089600  0.0044614  -2.008  0.04697 *  
## d(L(mfn_neg, 2))    -0.0062382  0.0039342  -1.586  0.11559    
## d(trend)                    NA         NA      NA       NA    
## d(L(trend, 1))              NA         NA      NA       NA    
## d(L(trend, 2))              NA         NA      NA       NA    
## month_1              0.0110494  0.0065907   1.677  0.09637 .  
## month_2              0.0151695  0.0073005   2.078  0.03997 *  
## month_3              0.0155581  0.0067551   2.303  0.02308 *  
## month_4              0.0138485  0.0076003   1.822  0.07106 .  
## month_5             -0.0001994  0.0086367  -0.023  0.98162    
## month_6             -0.0037108  0.0083474  -0.445  0.65749    
## month_7             -0.0062672  0.0087355  -0.717  0.47457    
## month_8             -0.0129618  0.0100075  -1.295  0.19786    
## month_9             -0.0054360  0.0086747  -0.627  0.53214    
## month_10            -0.0117755  0.0058816  -2.002  0.04765 *  
## month_11             0.0044759  0.0057177   0.783  0.43536    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005822 on 114 degrees of freedom
## Multiple R-squared:  0.9141, Adjusted R-squared:  0.8854 
## F-statistic: 31.91 on 38 and 114 DF,  p-value: < 2.2e-16

# ARDL modelinizin varsayımlarını kontrol etmek için:

# 1. Otokorelasyon testi
library(lmtest)
dwtest(modnardl$best_model) # Durbin-Watson testi

## 
##  Durbin-Watson test
## 
## data:  modnardl$best_model
## DW = 2.083, p-value = 0.7632
## alternative hypothesis: true autocorrelation is greater than 0

bgtest(modnardl$best_model, order = 12) # Breusch-Godfrey testi

## 
##  Breusch-Godfrey test for serial correlation of order up to 12
## 
## data:  modnardl$best_model
## LM test = 57.721, df = 12, p-value = 5.856e-08

# 2. Normallik testi
library(tseries)
jarque.bera.test(residuals(modnardl$best_model))

## 
##  Jarque Bera Test
## 
## data:  residuals(modnardl$best_model)
## X-squared = 27.36, df = 2, p-value = 1.145e-06

# 3. Değişen varyans testi
bptest(modnardl$best_model) # Breusch-Pagan testi

## 
##  studentized Breusch-Pagan test
## 
## data:  modnardl$best_model
## BP = 53.041, df = 38, p-value = 0.05333

library(sandwich)
coeftest(modnardl$best_model, vcov = vcovHC(modnardl$best_model, type = "HC0"))

## 
## t test of coefficients:
## 
##                     Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)       2.71718805  0.47233874  5.7526 7.504e-08 ***
## L(log_gdp, 1)     1.82119896  0.11388696 15.9913 < 2.2e-16 ***
## L(log_gdp, 2)    -1.43919839  0.18128736 -7.9388 1.570e-12 ***
## L(log_gdp, 3)     0.43105388  0.09415927  4.5779 1.205e-05 ***
## log_FA_pos        0.01227383  0.00715141  1.7163 0.0888270 .  
## L(log_FA_pos, 1)  0.00880427  0.00731879  1.2030 0.2314815    
## L(log_FA_pos, 2) -0.00564521  0.00479174 -1.1781 0.2412047    
## L(log_FA_pos, 3)  0.01371377  0.00377619  3.6316 0.0004234 ***
## log_FA_neg        0.01022953  0.00740486  1.3815 0.1698384    
## L(log_FA_neg, 1) -0.01356316  0.00838522 -1.6175 0.1085326    
## L(log_FA_neg, 2)  0.02355995  0.00774170  3.0433 0.0029058 ** 
## log_DA_pos       -0.00543295  0.00755462 -0.7192 0.4735163    
## L(log_DA_pos, 1)  0.01800298  0.01245929  1.4449 0.1512164    
## L(log_DA_pos, 2) -0.02114897  0.00762943 -2.7720 0.0065084 ** 
## log_DA_neg       -0.00225812  0.00696183 -0.3244 0.7462625    
## mdn_pos           0.00145586  0.00920867  0.1581 0.8746600    
## L(mdn_pos, 1)     0.02600702  0.01055336  2.4643 0.0152190 *  
## L(mdn_pos, 2)    -0.01324060  0.00675730 -1.9595 0.0525010 .  
## mdn_neg           0.01683713  0.01021631  1.6481 0.1020929    
## L(mdn_neg, 1)    -0.00732526  0.00746775 -0.9809 0.3287095    
## L(mdn_neg, 2)     0.01111264  0.00687272  1.6169 0.1086593    
## L(mdn_neg, 3)    -0.01295563  0.00624789 -2.0736 0.0403691 *  
## mfn_pos          -0.00098678  0.00408211 -0.2417 0.8094211    
## mfn_neg          -0.00093528  0.00536784 -0.1742 0.8619874    
## L(mfn_neg, 1)    -0.01192941  0.00569969 -2.0930 0.0385680 *  
## L(mfn_neg, 2)     0.00272183  0.00467046  0.5828 0.5611959    
## L(mfn_neg, 3)     0.00623817  0.00333080  1.8729 0.0636473 .  
## trend            -0.00111925  0.00068098 -1.6436 0.1030165    
## month_1           0.01104944  0.00519214  2.1281 0.0354819 *  
## month_2           0.01516946  0.00607013  2.4990 0.0138778 *  
## month_3           0.01555812  0.00568759  2.7355 0.0072264 ** 
## month_4           0.01384848  0.00559167  2.4766 0.0147312 *  
## month_5          -0.00019944  0.00738090 -0.0270 0.9784899    
## month_6          -0.00371085  0.00667005 -0.5563 0.5790658    
## month_7          -0.00626722  0.00814299 -0.7696 0.4431014    
## month_8          -0.01296182  0.00805205 -1.6098 0.1102171    
## month_9          -0.00543604  0.00714004 -0.7613 0.4480229    
## month_10         -0.01177545  0.00486260 -2.4216 0.0170274 *  
## month_11          0.00447588  0.00476032  0.9402 0.3490778    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 5. Model spesifikasyonu testi - RESET testi
resettest(modnardl$best_model)

## 
##  RESET test
## 
## data:  modnardl$best_model
## RESET = 0.50389, df1 = 2, df2 = 109, p-value = 0.6056

uecm_model <- uecm(modnardl$best_model)
# Uzun dönem ilişkiler
coef_uecm <- coef(uecm_model)

# log_FA için uzun dönem katsayılar
LR_log_FA_pos <- -coef_uecm["L(log_FA_pos, 1)"] / coef_uecm["L(log_gdp, 1)"]
LR_log_FA_neg <- -coef_uecm["L(log_FA_neg, 1)"] / coef_uecm["L(log_gdp, 1)"]

# log_DA için uzun dönem katsayılar
LR_log_DA_pos <- -coef_uecm["L(log_DA_pos, 1)"] / coef_uecm["L(log_gdp, 1)"]
LR_log_DA_neg <- -coef_uecm["log_DA_neg"] / coef_uecm["L(log_gdp, 1)"]

# mdn için uzun dönem katsayılar
LR_mdn_pos <- -coef_uecm["L(mdn_pos, 1)"] / coef_uecm["L(log_gdp, 1)"]
LR_mdn_neg <- -coef_uecm["L(mdn_neg, 1)"] / coef_uecm["L(log_gdp, 1)"]

# mfn için uzun dönem katsayılar
LR_mfn_pos <- -coef_uecm["mfn_pos"] / coef_uecm["L(log_gdp, 1)"]
LR_mfn_neg <- -coef_uecm["L(mfn_neg, 1)"] / coef_uecm["L(log_gdp, 1)"]

# Uzun dönem katsayıları gösterme
long_run_coefs <- data.frame(
  Positive = c(LR_log_FA_pos, LR_log_DA_pos, LR_mdn_pos, LR_mfn_pos),
  Negative = c(LR_log_FA_neg, LR_log_DA_neg, LR_mdn_neg, LR_mfn_neg),
  row.names = c("log_FA", "log_DA", "mdn", "mfn")
)
print(long_run_coefs)

##            Positive    Negative
## log_FA  0.155909859  0.10819366
## log_DA -0.045890039 -0.01207900
## mdn     0.076077126  0.04102195
## mfn    -0.005278446 -0.02088680

summary(recm(modnardl$best_model, case = 3))

## 
## Time series regression with "zooreg" data:
## Start = 4, End = 156
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0195675 -0.0027164 -0.0001845  0.0026585  0.0160229 
## 
## Coefficients: (3 not defined because of singularities)
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.7160688  0.3252831   8.350 1.19e-13 ***
## d(L(log_gdp, 1))     1.0081445  0.0599197  16.825  < 2e-16 ***
## d(L(log_gdp, 2))    -0.4310539  0.0623981  -6.908 2.33e-10 ***
## d(log_FA_pos)        0.0122738  0.0065111   1.885 0.061779 .  
## d(L(log_FA_pos, 1)) -0.0080686  0.0057426  -1.405 0.162531    
## d(L(log_FA_pos, 2)) -0.0137138  0.0029216  -4.694 7.03e-06 ***
## d(log_FA_neg)        0.0102295  0.0035169   2.909 0.004308 ** 
## d(L(log_FA_neg, 1)) -0.0235599  0.0064219  -3.669 0.000362 ***
## d(log_DA_pos)       -0.0054329  0.0082614  -0.658 0.512005    
## d(L(log_DA_pos, 1))  0.0211490  0.0074813   2.827 0.005487 ** 
## d(mdn_pos)           0.0014559  0.0066824   0.218 0.827896    
## d(L(mdn_pos, 1))     0.0132406  0.0070107   1.889 0.061295 .  
## d(mdn_neg)           0.0168371  0.0075035   2.244 0.026627 *  
## d(L(mdn_neg, 1))     0.0018430  0.0074042   0.249 0.803842    
## d(L(mdn_neg, 2))     0.0129556  0.0068643   1.887 0.061464 .  
## d(mfn_neg)          -0.0009353  0.0049963  -0.187 0.851815    
## d(L(mfn_neg, 1))    -0.0089600  0.0036099  -2.482 0.014410 *  
## d(L(mfn_neg, 2))    -0.0062382  0.0034289  -1.819 0.071296 .  
## d(trend)                    NA         NA      NA       NA    
## d(L(trend, 1))              NA         NA      NA       NA    
## d(L(trend, 2))              NA         NA      NA       NA    
## month_1              0.0110494  0.0045236   2.443 0.016003 *  
## month_2              0.0151695  0.0054036   2.807 0.005812 ** 
## month_3              0.0155581  0.0052065   2.988 0.003389 ** 
## month_4              0.0138485  0.0056078   2.470 0.014901 *  
## month_5             -0.0001994  0.0062377  -0.032 0.974545    
## month_6             -0.0037108  0.0060442  -0.614 0.540382    
## month_7             -0.0062672  0.0058143  -1.078 0.283192    
## month_8             -0.0129618  0.0066965  -1.936 0.055211 .  
## month_9             -0.0054360  0.0065262  -0.833 0.406481    
## month_10            -0.0117755  0.0041044  -2.869 0.004849 ** 
## month_11             0.0044759  0.0051517   0.869 0.386642    
## ect                 -0.1869455  0.0223676  -8.358 1.14e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005605 on 123 degrees of freedom
##   (0 observations deleted due to missingness)
## Multiple R-squared:  0.9141, Adjusted R-squared:  0.8938 
## F-statistic: 45.12 on 29 and 123 DF,  p-value: < 2.2e-16

library(kardl)

## 
## Attaching package: 'kardl'

## The following object is masked from 'package:nardl':
## 
##     cusum

kardl_m <- kardl(bdf, log_gdp ~  deterministic(month_1 + month_2 + month_3 + month_4 + month_5 + month_6 + month_7 + month_8 + month_9 + month_10 + month_11) + trend + asym(log_FA + log_DA + mfn + mdn), maxlag = 3, mode = "performance")

## Warning: The formula is not well constructed!
## Warning: The formula is not well constructed!

summary(kardl_m)

## ==============================================================
## KARDL Summary
## 
## 
## Call:
## lm(formula = as.formula(fmodel), data = data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0205210 -0.0027409  0.0000832  0.0028226  0.0171513 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.1567767  0.5890577   5.359 4.42e-07 ***
## L1.log_gdp      -0.2171924  0.0403845  -5.378 4.06e-07 ***
## L1.log_FA_POS    0.0318497  0.0076074   4.187 5.60e-05 ***
## L1.log_FA_NEG    0.0227102  0.0075620   3.003  0.00328 ** 
## L1.log_DA_POS   -0.0151745  0.0094678  -1.603  0.11176    
## L1.log_DA_NEG   -0.0057473  0.0098779  -0.582  0.56183    
## L1.mfn_POS       0.0027833  0.0087430   0.318  0.75080    
## L1.mfn_NEG      -0.0004796  0.0072554  -0.066  0.94741    
## L1.mdn_POS       0.0185070  0.0084895   2.180  0.03131 *  
## L1.mdn_NEG       0.0115092  0.0067807   1.697  0.09236 .  
## L1.d.log_gdp     0.9791312  0.0689010  14.211  < 2e-16 ***
## L2.d.log_gdp    -0.3686030  0.0888378  -4.149 6.46e-05 ***
## L0.d.log_FA_POS  0.0124750  0.0088873   1.404  0.16313    
## L1.d.log_FA_POS -0.0027550  0.0066337  -0.415  0.67870    
## L2.d.log_FA_POS -0.0093217  0.0040168  -2.321  0.02208 *  
## L0.d.log_FA_NEG  0.0166336  0.0085952   1.935  0.05544 .  
## L1.d.log_FA_NEG -0.0248954  0.0092249  -2.699  0.00802 ** 
## L2.d.log_FA_NEG -0.0094804  0.0059877  -1.583  0.11612    
## L0.d.log_DA_POS -0.0058274  0.0097184  -0.600  0.54995    
## L1.d.log_DA_POS  0.0279364  0.0087416   3.196  0.00180 ** 
## L0.d.log_DA_NEG -0.0064787  0.0112662  -0.575  0.56639    
## L1.d.log_DA_NEG -0.0051765  0.0110069  -0.470  0.63904    
## L0.d.mfn_POS    -0.0007944  0.0044859  -0.177  0.85975    
## L0.d.mfn_NEG     0.0077181  0.0088521   0.872  0.38510    
## L0.d.mdn_POS     0.0043718  0.0090175   0.485  0.62874    
## L0.d.mdn_NEG     0.0129567  0.0095383   1.358  0.17702    
## month_1          0.0056796  0.0063720   0.891  0.37463    
## month_2          0.0167389  0.0078321   2.137  0.03472 *  
## month_3          0.0183885  0.0070413   2.612  0.01023 *  
## month_4          0.0170289  0.0080577   2.113  0.03675 *  
## month_5          0.0022002  0.0084043   0.262  0.79395    
## month_6         -0.0011582  0.0079511  -0.146  0.88444    
## month_7         -0.0049547  0.0087665  -0.565  0.57306    
## month_8         -0.0110262  0.0103042  -1.070  0.28685    
## month_9         -0.0055201  0.0087065  -0.634  0.52734    
## month_10        -0.0086106  0.0059789  -1.440  0.15256    
## month_11         0.0052230  0.0052235   1.000  0.31948    
## trend           -0.0008776  0.0006340  -1.384  0.16897    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005978 on 114 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.9024, Adjusted R-squared:  0.8707 
## F-statistic: 28.49 on 37 and 114 DF,  p-value: < 2.2e-16
## 
## ===============================================================

MyCoint<- aardl(kardl_m)
# critical Values are
MyCoint$criticalValues

##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   7.00   1.72   3.17   2.02   3.62   2.31   4.04   2.68   4.54

mySummary<-summary(MyCoint)
 mySummary

## Method: AARDL
## H0: Coef(L1.log_FA_POS) = Coef(L1.log_FA_NEG) = Coef(L1.log_DA_POS) = Coef(L1.log_DA_NEG) = Coef(L1.mfn_POS) = Coef(L1.mfn_NEG) = Coef(L1.mdn_POS) = Coef(L1.mdn_NEG) = 0
## H1: Coef(L1.log_FA_POS) ≠ Coef(L1.log_FA_NEG) ≠ Coef(L1.log_DA_POS) ≠ Coef(L1.log_DA_NEG) ≠ Coef(L1.mfn_POS) ≠ Coef(L1.mfn_NEG) ≠ Coef(L1.mdn_POS) ≠ Coef(L1.mdn_NEG)≠ 0
## AARDL test.
## The calculated AARDL Wald F is 4.84765886082129 and  k=7
##  The calculated Wald F for PSS is 5.60573698142598
##  The calculated PSS t is -5.37811300485268
##   
## The results: There is Cointegration 
## Significance level: 0.025 
## ___________________________________________________
## Critique values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   7.00   1.72   3.17   2.02   3.62   2.31   4.04   2.68   4.54

 mySummary$H0

## [1] "H0: Coef(L1.log_FA_POS) = Coef(L1.log_FA_NEG) = Coef(L1.log_DA_POS) = Coef(L1.log_DA_NEG) = Coef(L1.mfn_POS) = Coef(L1.mfn_NEG) = Coef(L1.mdn_POS) = Coef(L1.mdn_NEG) = 0"

# Summary of ARCH test
summary(archtest(kardl_m$finalModel$model$residuals))

## ARCH test
## 
## F = 3.236653, p-value = 0.07403146, df1 = 1, df2 = 149
## 
## 
## Call:
## lm(formula = as.formula(paste0("resid~", paste("resid", 1:q, 
##     sep = "", collapse = "+"))), data = ndata)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -5.308e-05 -2.234e-05 -1.723e-05 -8.700e-07  3.843e-04 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.283e-05  4.918e-06   4.643 7.48e-06 ***
## resid1      1.459e-01  8.110e-02   1.799    0.074 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.413e-05 on 149 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.02126,    Adjusted R-squared:  0.01469 
## F-statistic: 3.237 on 1 and 149 DF,  p-value: 0.07403

ast<- asymmetrytest(kardl_m)
 ast

## Asymmetries in the long run
## 
##                F          P
## log_FA 4.4638621 0.03680132
## log_DA 2.1784140 0.14271546
## mfn    0.7346251 0.39318631
## mdn    1.0271164 0.31298308
## ________________________________
## Asymmetries in the short run
## 
##                F         P
## log_FA 0.8177091 0.3677578
## log_DA 2.5624004 0.1121999
## mfn    0.7935483 0.3749056
## mdn    0.4102343 0.5231360

 ast$Lhypotheses

## $H0
## [1] "- Coef(L1.log_FA_POS)/Coef(L1.log_gdp) = - Coef(L1.log_FA_NEG)/Coef(L1.log_gdp)"
## [2] "- Coef(L1.log_DA_POS)/Coef(L1.log_gdp) = - Coef(L1.log_DA_NEG)/Coef(L1.log_gdp)"
## [3] "- Coef(L1.mfn_POS)/Coef(L1.log_gdp) = - Coef(L1.mfn_NEG)/Coef(L1.log_gdp)"      
## [4] "- Coef(L1.mdn_POS)/Coef(L1.log_gdp) = - Coef(L1.mdn_NEG)/Coef(L1.log_gdp)"      
## 
## $H1
## [1] "- Coef(L1.log_FA_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.log_FA_NEG)/Coef(L1.log_gdp)"
## [2] "- Coef(L1.log_DA_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.log_DA_NEG)/Coef(L1.log_gdp)"
## [3] "- Coef(L1.mfn_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.mfn_NEG)/Coef(L1.log_gdp)"      
## [4] "- Coef(L1.mdn_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.mdn_NEG)/Coef(L1.log_gdp)"

 summary(ast)

## Asymmetries in the long run
## H0: - Coef(L1.log_FA_POS)/Coef(L1.log_gdp) = - Coef(L1.log_FA_NEG)/Coef(L1.log_gdp)
## H1: - Coef(L1.log_FA_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.log_FA_NEG)/Coef(L1.log_gdp)
## 
## H0: - Coef(L1.log_DA_POS)/Coef(L1.log_gdp) = - Coef(L1.log_DA_NEG)/Coef(L1.log_gdp)
## H1: - Coef(L1.log_DA_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.log_DA_NEG)/Coef(L1.log_gdp)
## 
## H0: - Coef(L1.mfn_POS)/Coef(L1.log_gdp) = - Coef(L1.mfn_NEG)/Coef(L1.log_gdp)
## H1: - Coef(L1.mfn_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.mfn_NEG)/Coef(L1.log_gdp)
## 
## H0: - Coef(L1.mdn_POS)/Coef(L1.log_gdp) = - Coef(L1.mdn_NEG)/Coef(L1.log_gdp)
## H1: - Coef(L1.mdn_POS)/Coef(L1.log_gdp) ≠ - Coef(L1.mdn_NEG)/Coef(L1.log_gdp)
## 
##                F          P df1 df2
## log_FA 4.4638621 0.03680132   1 114
## log_DA 2.1784140 0.14271546   1 114
## mfn    0.7346251 0.39318631   1 114
## mdn    1.0271164 0.31298308   1 114
## _____________________________
## 
## Asymmetries in the long run
## H0: Coef(L0.d.log_FA_POS) + Coef(L1.d.log_FA_POS) + Coef(L2.d.log_FA_POS) = Coef(L0.d.log_FA_NEG) + Coef(L1.d.log_FA_NEG) + Coef(L2.d.log_FA_NEG)
## H1: Coef(L0.d.log_FA_POS) + Coef(L1.d.log_FA_POS) + Coef(L2.d.log_FA_POS) ≠ Coef(L0.d.log_FA_NEG) + Coef(L1.d.log_FA_NEG) + Coef(L2.d.log_FA_NEG)
## 
## H0: Coef(L0.d.log_DA_POS) + Coef(L1.d.log_DA_POS) = Coef(L0.d.log_DA_NEG) + Coef(L1.d.log_DA_NEG)
## H1: Coef(L0.d.log_DA_POS) + Coef(L1.d.log_DA_POS) ≠ Coef(L0.d.log_DA_NEG) + Coef(L1.d.log_DA_NEG)
## 
## H0: Coef(L0.d.mfn_POS) = Coef(L0.d.mfn_NEG)
## H1: Coef(L0.d.mfn_POS) ≠ Coef(L0.d.mfn_NEG)
## 
## H0: Coef(L0.d.mdn_POS) = Coef(L0.d.mdn_NEG)
## H1: Coef(L0.d.mdn_POS) ≠ Coef(L0.d.mdn_NEG)
## 
##                F         P DF    Sum.of.Sq ResDF1 ResDF2        RSS1
## log_FA 0.8177091 0.3677578  1 2.921967e-05    115    114 0.004102848
## log_DA 2.5624004 0.1121999  1 9.156374e-05    115    114 0.004165192
## mfn    0.7935483 0.3749056  1 2.835632e-05    115    114 0.004101984
## mdn    0.4102343 0.5231360  1 1.465914e-05    115    114 0.004088287
##               RSS2
## log_FA 0.004073628
## log_DA 0.004073628
## mfn    0.004073628
## mdn    0.004073628

MyCoint<- ecmtest(kardl_m)

## Warning: Attention: ECMtest function is deprecated.

cat(paste0("The t statistics=",MyCoint$result$t," where k=",MyCoint$result$k,". \nWe found '",MyCoint$result$Cont, "' at ",names(MyCoint$result$star),"."))

## The t statistics= where k=. 
## We found '' at .

# critical Values are
MyCoint$criticalValues

##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   7.00  -3.13  -4.53  -3.41  -4.85  -3.65  -5.14  -3.96  -5.49

mySummary<-summary(MyCoint)
 mySummary

## Method: ECM
## H0: Coef(L1.log_gdp) = 0
## H1: Coef(L1.log_gdp)≠ 0
## Attention: this function is deprecated.
## The coeffient of ECM is -0.206672540946615 , t=-5.70965218593432 k=8  
## The results: There is Cointegration 
## Significance level: 0.01 
## ___________________________________________________
## Critique values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   7.00  -3.13  -4.53  -3.41  -4.85  -3.65  -5.14  -3.96  -5.49

m<-mplier(kardl_m,40)
 head(m$mpsi)

##      horizon log_FA_POS  log_FA_NEG   log_FA_dif   log_DA_POS  log_DA_NEG
## [1,]       0 0.01247501 -0.01663362 -0.004158607 -0.005827351 0.006478672
## [2,]       1 0.05107496 -0.02712214  0.023952827  0.002494386 0.022338824
## [3,]       2 0.09570597 -0.03859958  0.057106392 -0.002925862 0.036375388
## [4,]       3 0.13624064 -0.06029802  0.075942613 -0.025839475 0.042119745
## [5,]       4 0.16173757 -0.08692692  0.074810649 -0.055839390 0.039169492
## [6,]       5 0.16848277 -0.10883229  0.059650481 -0.079813865 0.031403381
##         log_DA_dif       mfn_POS      mfn_NEG      mfn_dif    mdn_POS
## [1,]  0.0006513208 -0.0007944467 -0.007718104 -0.008512551 0.00437182
## [2,]  0.0248332097  0.0013835356 -0.013119204 -0.011735669 0.02620985
## [3,]  0.0334495262  0.0062917104 -0.012233660 -0.005941949 0.05879508
## [4,]  0.0162802696  0.0117114374 -0.006239053  0.005472384 0.08838786
## [5,] -0.0166698979  0.0154485607  0.001138736  0.016587297 0.10466187
## [6,] -0.0484104839  0.0165379601  0.006385229  0.022923189 0.10546349
##          mdn_NEG      mdn_dif
## [1,] -0.01295668 -0.008584865
## [2,] -0.03433807 -0.008128215
## [3,] -0.05454859  0.004246491
## [4,] -0.06611775  0.022270115
## [5,] -0.06714473  0.037517147
## [6,] -0.06081171  0.044651784

 head(m$omega)

## [1]  1.761939 -1.347734  0.368603

 head(m$Lambda)

##        log_FA_POS   log_FA_NEG   log_DA_POS   log_DA_NEG       mfn_POS
## [1,]  0.012475013  0.016633620 -0.005827351 -0.006478672 -0.0007944467
## [2,]  0.016619743 -0.018818902  0.018589173 -0.004445128  0.0035777490
## [3,] -0.006566749  0.015415011 -0.027936360  0.005176526  0.0000000000
## [4,]  0.009321715  0.009480432  0.000000000  0.000000000  0.0000000000
## [5,]  0.000000000  0.000000000  0.000000000  0.000000000  0.0000000000
## [6,]  0.000000000  0.000000000  0.000000000  0.000000000  0.0000000000
##           mfn_NEG    mdn_POS      mdn_NEG
## [1,]  0.007718104 0.00437182  0.012956685
## [2,] -0.008197727 0.01413515 -0.001447502
## [3,]  0.000000000 0.00000000  0.000000000
## [4,]  0.000000000 0.00000000  0.000000000
## [5,]  0.000000000 0.00000000  0.000000000
## [6,]  0.000000000 0.00000000  0.000000000

 mydata<-as.data.frame( m$mpsi)
 ggplot(mydata, aes(x = 1:nrow(mydata))) +
   geom_line(aes(y = log_FA_POS  ), color = "blue", linetype = "solid") +
   geom_line(aes(y = log_FA_NEG   ), color = "red", linetype = "dashed") +
   geom_line(aes(y = log_FA_dif   ), color = "black", linetype = "dotted") +
   labs(x = "Time", y = "Value", title = "Dynamic multipliers with ggplot2") +
   theme_minimal()

ggplot(mydata, aes(x = 1:nrow(mydata))) +
   geom_line(aes(y = mfn_POS  ), color = "blue", linetype = "solid") +
   geom_line(aes(y = mfn_NEG    ), color = "red", linetype = "dashed") +
   geom_line(aes(y = mfn_dif    ), color = "black", linetype = "dotted") +
   labs(x = "Time", y = "Value", title = "Dynamic multipliers with ggplot2") +
   theme_minimal()

MyCoint<-narayan(kardl_m)

## Warning: Warning: Narayan F test is designated for small periods. Your data's length is 155
## You can use pssf test for this purpose. However, the calculation has been made here.

MyCoint

## Narayan F test.
## The F statistics = 5.60573698142598 where k = 4. 
## The proboblity is = 0.0000021510

mySummary<-summary(MyCoint)
mySummary

## Method: Narayan
## H0: Coef(L1.log_gdp) = Coef(L1.log_FA_POS) = Coef(L1.log_FA_NEG) = Coef(L1.log_DA_POS) = Coef(L1.log_DA_NEG) = Coef(L1.mfn_POS) = Coef(L1.mfn_NEG) = Coef(L1.mdn_POS) = Coef(L1.mdn_NEG) = 0
## H1: Coef(L1.log_gdp) ≠ Coef(L1.log_FA_POS) ≠ Coef(L1.log_FA_NEG) ≠ Coef(L1.log_DA_POS) ≠ Coef(L1.log_DA_NEG) ≠ Coef(L1.mfn_POS) ≠ Coef(L1.mfn_NEG) ≠ Coef(L1.mdn_POS) ≠ Coef(L1.mdn_NEG)≠ 0
## The F statistics = 5.60573698142598 where k = 4.
## The proboblity is = 0.0000021510  
## The results: There is Cointegration 
## Significance level: 0.05 
## ___________________________________________________
## Critique values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##  4.000  3.160  4.230  3.678  4.840     NA     NA  4.890  6.164

MyCoint<-pssf(kardl_m)
summary(MyCoint)

## Method: PesaranF
## H0: Coef(L1.log_gdp) = Coef(L1.log_FA_POS) = Coef(L1.log_FA_NEG) = Coef(L1.log_DA_POS) = Coef(L1.log_DA_NEG) = Coef(L1.mfn_POS) = Coef(L1.mfn_NEG) = Coef(L1.mdn_POS) = Coef(L1.mdn_NEG) = 0
## H1: Coef(L1.log_gdp) ≠ Coef(L1.log_FA_POS) ≠ Coef(L1.log_FA_NEG) ≠ Coef(L1.log_DA_POS) ≠ Coef(L1.log_DA_NEG) ≠ Coef(L1.mfn_POS) ≠ Coef(L1.mfn_NEG) ≠ Coef(L1.mdn_POS) ≠ Coef(L1.mdn_NEG)≠ 0
## The F statistics = 5.60573698142598 where k = 4. 
## The proboblity is = 0.0000021510  
## The results: There is Cointegration 
## Significance level: 0.025 
## ___________________________________________________
## Critique values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   4.00   3.03   4.06   3.47   4.57   3.89   5.07   4.40   5.72

MyCoint<-psst(kardl_m)
summary(MyCoint)

## Method: Pesarant
## H0: Coef(L1.log_gdp) = 0
## H1: Coef(L1.log_gdp)≠ 0
## The t statistics = -5.37811300485268 where k = 8.  
## The results: There is Cointegration 
## Significance level: 0.025 
## ___________________________________________________
## Critique values are as follows:
##      k  0.10L  0.10U  0.05L  0.05U 0.025L 0.025U  0.01L  0.01U 
##   7.00  -3.13  -4.53  -3.41  -4.85  -3.65  -5.14  -3.96  -5.49

Otoregresif Dağıtılmış Gecikme (ARDL) Modelleri: Teori

EViews Blog

2025-03-08

Giriş

Modelin Belirlenmesi (Specification)

İlk Temsil: Zamana Yayılmış Dinamik Regresyon (Intertemporal Dynamics Regression)

İkinci Temsil: Uzun Dönem Dinamiklerinin Regresyon Sonrası Türetilmesi

Üçüncü Temsil: Koşullu Hata Düzeltme Formu ve Sınır Testi

Case 1: (Sabit Terim ve Trend Yok): \(a_0 = a_1 = 0\), yani, \((\boldsymbol{\mu} = \boldsymbol{\gamma} = \mathbf{0})\)

DGP:

Regresyon:

Case 2: (Kısıtlı Sabit Terim ve Trend Yok): \(a_0 = \psi(1)\mu_y + \sum_{j=1}^k \beta_j(1)\mu_{x_j}\) ve \(a_1 = 0\) öyle ki \(\boldsymbol{\gamma} = \mathbf{0}\)

DGP:

Regresyon:

Durum 3: (Kısıtlanmamış Sabit ve Trend Yok): \(a_0 \neq 0\) ve \(a_1 = 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(\boldsymbol{\gamma} = \mathbf{0}\)

DGP:

Regresyon:

Durum 4: (Kısıtlanmamış Sabit ve Kısıtlı Trend): \(a_0 \neq 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(a_1 = \psi(1)\gamma_y + \sum_{j=1}^{k} \beta_j(1)\gamma_{x_j}\)

DGP:

Regresyon:

Durum 5: (Kısıtlanmamış Sabit ve Kısıtlanmamış Trend): \(a_0 \neq 0\) ve \(a_1 \neq 0\) öyle ki \(\boldsymbol{\mu} \neq \mathbf{0}\) ve \(\boldsymbol{\gamma} \neq \mathbf{0}\)

DGP:

Regresyon:

ARDL Model Uygulama Prosedürü

Faiz Oranı Termin Yapısı (TSIR) ve Eşbütünleşme

Deterministik Belirlemeler

ARDL Gecikme Yapısının Belirlenmesi

Model 1: Eşbütünleşme yok (No Cointegrating Relationship)

NARDL: A Study of Bosnia’s Tourism Sector

Veri ve Motivasyon

Asimetri Tahmini