Mô hình tự hồi quy véc-tơ là một công cụ mạnh mẽ và phổ biến để phân tích động thái của nhiều chuỗi thời gian cùng lúc. Mô hình này biểu diễn một véc-tơ gồm các biến quan sát dưới dạng một hàm số theo các trễ của chính chúng. Trong bài viết này, chúng ta sẽ tìm hiểu cách mô phỏng dữ liệu từ một quy trình tự hồi quy véc-tơ bậc hai với hai biến số, tiến hành ước lượng tham số bằng Stata và thực hiện suy diễn thống kê thông qua các hàm phản ứng đẩy.
Mô phỏng dữ liệu
Hãy bắt đầu bằng việc mô phỏng một quy trình tự hồi quy véc-tơ bậc hai với hai biến số. Phương trình tổng quát có dạng hệ phương trình tuyến tính, trong đó các biến hiện tại phụ thuộc vào giá trị trễ bậc một và bậc hai của chính chúng và biến còn lại, cộng với một véc-tơ nhiễu trắng tuân theo phân phối chuẩn đa biến với trung bình bằng không và ma trận hiệp phương sai cho trước.
Chúng ta thiết lập kích thước mẫu ban đầu là 1100 quan sát, sau đó tạo biến thời gian và các biến chứa chuỗi quan sát cũng như các biến nhiễu.
1clear all
2set seed 2016
3local T = 1100
4set obs `T'
5gen time = _n
6tsset time
7generate y1 = .
8generate y2 = .
9generate eps1 = .
10generate eps2 = .Trong đoạn mã trên, chúng ta đặt hạt giống ngẫu nhiên để đảm bảo tính tái lập, thiết lập kích thước mẫu và định dạng biến thời gian chuỗi số. Tiếp theo, các biến trống được khởi tạo để sẵn sàng nhận dữ liệu mô phỏng.
Thiết lập giá trị tham số
Chúng ta chọn các giá trị tham số cụ thể cho mô hình tự hồi quy véc-tơ bậc hai. Hệ số chặn của chuỗi thứ nhất là 0.1, chuỗi thứ hai là 0.4. Các ma trận hệ số trễ bậc một và bậc hai, cùng ma trận hiệp phương sai của sai số được định nghĩa trong môi trường Mata của Stata như sau:
1mata:
2mu = (0.1\0.4)
3A1 = (0.6,-0.3\0.4,0.2)
4A2 = (0.2,0.3\-0.1,0.1)
5Sigma = (1,0.5\0.5,1)end
Trước khi tạo dữ liệu, việc kiểm tra xem các tham số đã chọn có thỏa mãn điều kiện ổn định của mô hình tự hồi quy véc-tơ hay không là rất quan trọng. Một quy trình tự hồi quy véc-tơ bậc hai ổn định nếu tất cả các trị riêng của ma trận đồng hành F đều có mô-đun nhỏ hơn 1. Chúng ta tính toán các trị riêng này bằng đoạn mã sau:
1mata:
2K = p = 2
3F = J(K*p,K*p,0)
4F[1..2,1..2] = A1
5F[1..2,3..4] = A2
6F[3..4,1..2] = I(K)
7X = L = .
8eigensystem(F,X,L)
9L'end
Kết quả tính toán cho thấy tất cả các trị riêng đều nằm trong đường tròn đơn vị, tức là có giá trị tuyệt đối nhỏ hơn 1. Điều này đảm bảo quy trình mô phỏng sẽ ổn định và không bị bùng nổ theo thời gian.
Tạo sai số ngẫu nhiên từ phân phối chuẩn đa biến
Bước tiếp theo là rút các giá trị ngẫu nhiên từ phân phối chuẩn đa biến với ma trận hiệp phương sai Sigma để gán vào các biến sai số trong Stata.
1mata:
2T = strtoreal(st_local("T"))
3u = rnormal(T,2,0,1)*cholesky(Sigma)
4epsmat = .
5st_view(epsmat,.,"eps1 eps2")
6epsmat[1..T,.] = uend
Chúng ta chuyển kích thước mẫu từ biến cục bộ của Stata sang Mata, sau đó sử dụng phép phân tích Cholesky để chuyển đổi các biến ngẫu nhiên chuẩn tắc độc lập thành các biến ngẫu nhiên có cấu trúc tương quan mong muốn. Hàm st_view giúp đẩy trực tiếp các giá trị này vào bộ dữ liệu Stata hiện hành.
Tạo chuỗi quan sát thực tế
Để các quan sát đầu tiên có cùng cấu trúc tương quan với phần còn lại của mẫu, chúng ta áp dụng phương pháp của Lutkepohl. Các giá trị khởi đầu được tạo ra dựa trên kỳ vọng vô điều kiện và ma trận hiệp phương sai dài hạn của chuỗi.
1mata:
2Sigma_e = J(K*p,K*p,0)
3Sigma_e[1..K,1..K] = Sigma
4Sigma_y = luinv(I((K*p)^2)-F#F)*vec(Sigma_e)
5Sigma_y = rowshape(Sigma_y,K*p)'
6theta = luinv(I(K)-A1-A2)*mu
7Q = cholesky(Sigma_y)*rnormal(K*p,1,0,1)
8data = .
9st_view(data,.,"y1 y2")
10data[1..p,.] = ((Q[3..4],Q[1..2]):+mu)'end
Sau khi đã có hai quan sát khởi đầu, chúng ta dùng vòng lặp trong Stata để tạo ra toàn bộ chuỗi thời gian theo đúng phương trình sai phân của mô hình tự hồi quy véc-tơ bậc hai. Chúng ta cũng loại bỏ 100 quan sát đầu tiên để xóa bỏ ảnh hưởng của các giá trị khởi đầu.
1forvalues i=3/`T' {
2 qui {
3replace y1 = 0.1 + 0.6*l.y1 - 0.3*l.y2 + 0.2*l2.y1 + 0.3*l2.y2 + eps1 in `i'
4replace y2 = 0.4 + 0.4*l.y1 + 0.2*l.y2 - 0.1*l2.y1 + 0.1*l2.y2 + eps2 in `i'
5 }
6}
7drop in 1/100Ước lượng mô hình
Sau khi đã có dữ liệu mô phỏng, chúng ta tiến hành ước lượng các tham số bằng lệnh var trong Stata để kiểm tra xem mô hình ước lượng có tiệm cận gần với các tham số gốc ban đầu hay không.
1var y1 y2Kết quả ước lượng cho thấy các hệ số thu được rất gần với các giá trị thực tế mà chúng ta đã thiết lập lúc đầu. Các sai số tiêu chuẩn nhỏ và các hệ số đều có ý nghĩa thống kê cao, chứng minh rằng phương pháp mô phỏng và ước lượng hoạt động cực kỳ chính xác.
Suy diễn thống kê qua hàm phản ứng đẩy
Hàm phản ứng đẩy là công cụ đắc lực để phân tích tác động của một cú sốc ngoại sinh từ sai số của một phương trình lên các biến nội sinh trong hệ thống theo thời gian. Để tính toán hàm phản ứng đẩy, mô hình tự hồi quy véc-tơ bậc hai được chuyển đổi sang dạng trung bình trượt vô hạn.
Các phản ứng trong tương lai của chuỗi do một cú sốc đơn vị gây ra chính là các hệ số của ma trận trung bình trượt này. Chúng ta có thể tạo và lưu trữ kết quả tính toán hàm phản ứng đẩy bằng lệnh sau:
1irf create firstirf, set(myirf)Sau khi tạo, chúng ta có thể hiển thị bảng kết quả phản ứng của biến y1 khi có cú sốc xuất hiện từ chính phương trình của nó bằng cách gõ:
1irf table irf, impulse(y1) response(y1) nociBảng kết quả hiển thị mức độ ảnh hưởng giảm dần theo thời gian, quay trở về trạng thái cân bằng dài hạn. Để trực quan hóa động thái này cùng với khoảng tin cậy 95 phần trăm, chúng ta vẽ biểu đồ phản ứng đẩy.
1irf graph irf, impulse(y1) response(y1)
Biểu đồ cho thấy khi có một cú sốc đơn vị xảy ra, biến y1 tăng ngay lập tức một đơn vị tại thời điểm ban đầu, sau đó giảm dần qua các bước thời gian tiếp theo về mức cân bằng ban đầu.
Hàm phản ứng đẩy trực giao
Trong thực tế, các sai số giữa các phương trình thường có sự tương quan mạnh mẽ với nhau, nghĩa là chúng ta khó có thể giả định một cú sốc xảy ra ở phương trình này mà không ảnh hưởng lập tức đến phương trình khác. Hàm phản ứng đẩy trực giao giải quyết vấn đề này bằng cách sử dụng phép phân tích Cholesky để phân rã ma trận hiệp phương sai thành một ma trận tam giác dưới, từ đó cô lập các cú sốc độc lập với nhau.
Chúng ta có thể kiểm tra ma trận hiệp phương sai ước lượng được để thấy rõ mối tương quan này:
1matrix Sigma_hat = e(Sigma)
2matrix P_hat = cholesky(Sigma_hat)
3matrix list P_hatBằng cách nhân ma trận hệ số trung bình trượt với ma trận tam giác dưới này, chúng ta thu được hàm phản ứng đẩy trực giao. Chúng ta hiển thị bảng kết quả của các hàm phản ứng trực giao bằng lệnh sau:
1irf table oirf, nociCuối cùng, chúng ta vẽ biểu đồ hàm phản ứng đẩy trực giao để thấy rõ tác động của cú sốc đã được loại bỏ sự đồng biến giữa các sai số.
1irf graph oirf, impulse(y1) response(y1)
Biểu đồ này giúp các nhà nghiên cứu đưa ra các suy diễn chính xác hơn về mặt kinh tế học hoặc các hiện tượng thực tế, khi các cú sốc không còn bị nhiễu bởi các yếu tố tác động đồng thời không quan sát được.
Kết luận
Bài viết đã hướng dẫn chi tiết cách thực hiện toàn bộ quy trình từ mô phỏng dữ liệu của một hệ tự hồi quy véc-tơ ổn định, ước lượng các tham số bằng công cụ chuyên dụng trong Stata, cho đến cách tính toán và biểu diễn đồ thị các hàm phản ứng đẩy thông thường cũng như trực giao. Việc hiểu rõ sự khác biệt giữa hai loại hàm phản ứng này là chìa khóa để đưa ra các phân tích định lượng chính xác trong nghiên cứu kinh tế lượng chuỗi thời gian.
✨ Điểm mấu chốt khi phân tích mô hình tự hồi quy véc-tơ là việc lựa chọn giữa hàm phản ứng đẩy thông thường và hàm phản ứng đẩy trực giao. Nếu ma trận hiệp phương sai của sai số có mối tương quan đồng thời lớn, hàm trực giao phân rã Cholesky sẽ là lựa chọn bắt buộc để tránh các suy diễn sai lệch về tác động của chính sách hoặc các cú sốc thị trường.
Nếu thứ tự sắp xếp các biến trong mô hình tự hồi quy véc-tơ thay đổi, kết quả của hàm phản ứng đẩy trực giao sẽ bị ảnh hưởng như thế nào? Bạn hãy thử hoán đổi vị trí của hai biến y1 và y2 trong câu lệnh ước lượng Cholesky trên Stata và so sánh đồ thị nhận được để tự tìm câu trả lời.
