Tiểu Luận Các mô hình ngẫu nhiên

Lời cảm ơn

Kính gửi Thầy Nguyễn Tuấn Duy Em Phạm Xuân Hoan xin gửi lời cảm ơn chân thành và sâu sắc nhất đến Thầy đã tận tình hướng dẫn em hoàn thành bài tiểu luận Phân tích sự phụ thuộc của TTCK Việt Nam đối với sự biến động của TTCK Nhật Bản với phương pháp copula có điều kiện . Sự tận tâm, nhiệt tình và những lời khuyên quý báu của Thầy đã giúp em rất nhiều trong quá trình nghiên cứu và hoàn thiện bài viết. Nhờ sự hướng dẫn của Thầy, em đã có thể hiểu rõ hơn về đề tài, trau dồi kiến thức và kỹ năng nghiên cứu, cũng như trình bày bài viết một cách khoa học và logic. Em xin trân trọng cảm ơn những góp ý của Thầy đã giúp em hoàn thiện bài tiểu luận một cách tốt nhất.

Lời cam kết

Em xin cam đoan số liệu và kết quả nghiên cứu trong báo cáo thực hành nghề nghiệp này là trung thực và các thông tin trích dẫn trong báo cáo đã được chỉ rõ nguồn gốc rõ ràng và được phép công bố.

Chương 1: Giới thiệu

1.1 Đặt vấn đề

Việc nghiên cứu cấu trúc phụ thuộc giữa thị trường chứng khoán Việt Nam và Nhật Bản thông qua phương pháp Copula có điều kiện là một hướng đi nghiên cứu hết sức cấp thiết trong bối cảnh hội nhập kinh tế toàn cầu ngày càng sâu rộng. Phương pháp Copula, với khả năng mô hình hóa linh hoạt các mối quan hệ phi tuyến tính và các cấu trúc phụ thuộc phức tạp, đặc biệt phù hợp để nghiên cứu mối quan hệ giữa hai thị trường chứng khoán có tính chất khác nhau như Việt Nam và Nhật Bản. Bằng cách áp dụng phương pháp này, chúng ta có thể khám phá những thay đổi trong cấu trúc phụ thuộc theo thời gian, chẳng hạn như sự gia tăng hoặc giảm đi của mối liên kết giữa hai thị trường trong các giai đoạn khủng hoảng hoặc phục hồi kinh tế.

Kết quả nghiên cứu không chỉ cung cấp cái nhìn sâu sắc về mối quan hệ giữa hai nền kinh tế mà còn có ý nghĩa thực tiễn to lớn. Đối với các nhà đầu tư, việc hiểu rõ cấu trúc phụ thuộc này giúp họ xây dựng danh mục đầu tư đa dạng hóa rủi ro hiệu quả hơn, từ đó nâng cao lợi nhuận và giảm thiểu tổn thất. Đối với các nhà hoạch định chính sách, nghiên cứu này cung cấp cơ sở khoa học để đưa ra các quyết sách điều tiết thị trường phù hợp, đảm bảo sự ổn định và phát triển bền vững của thị trường chứng khoán Việt Nam. Ngoài ra, việc áp dụng phương pháp Copula có điều kiện còn mở ra nhiều hướng nghiên cứu mới, như việc so sánh hiệu quả của các hàm Copula khác nhau, hoặc việc kết hợp phương pháp Copula với các mô hình biến động (GARCH) để nắm bắt sự biến động của thị trường một cách toàn diện hơn.

Với đà tăng trưởng xuất khẩu ấn tượng, gần gấp 8 lần kể từ khi gia nhập WTO, Việt Nam đã nhanh chóng vươn mình trở thành một trong những cường quốc thương mại hàng đầu thế giới. Song song đó, dòng vốn đầu tư trực tiếp nước ngoài (FDI), đặc biệt từ Nhật Bản, liên tục đổ vào Việt Nam, khẳng định vị thế hấp dẫn của thị trường nước ta. Tuy nhiên, sự hội nhập sâu rộng vào nền kinh tế toàn cầu cũng đồng nghĩa với việc thị trường Việt Nam dễ bị tác động bởi những biến động trên thị trường quốc tế. Nhận thức rõ tầm quan trọng của việc nắm bắt mối liên hệ giữa thị trường Việt Nam và các thị trường đối tác lớn, nghiên cứu này tập trung vào việc đánh giá cấu trúc phụ thuộc giữa thị trường chứng khoán Việt Nam và Nhật Bản. Việc chọn Nhật Bản là đối tượng nghiên cứu dựa trên thực tế là nước này là một trong những nhà đầu tư FDI lớn nhất vào Việt Nam. Bằng cách áp dụng phương pháp Copula có điều kiện, nghiên cứu sẽ giúp chúng ta hiểu rõ hơn về mức độ và tính chất của mối quan hệ này, từ đó hỗ trợ các nhà đầu tư, nhà hoạch định chính sách và các tổ chức tài chính trong việc ra quyết định, quản lý rủi ro và xây dựng các mô hình dự báo chính xác hơn trong một môi trường thị trường đầy biến động. Đó là lý do bài tiểu luận này nghiên cứu về “Đánh giá cấu trúc phụ thuộc của TTCK Việt Nam so với sự biến động của TTCK Nhật với phương pháp Copula có điều kiện”

1.2 Mục tiêu nghiên cứu

1.2.1 Mục tiêu tổng quát

Nghiên cứu mối tương quan và mức độ ảnh hưởng của TTCK Việt Nam đối với biến động của TTCK Nhật Bản trong giai đoạn 04/2018 - 7/2024 để đưa ra các đề xuất những kiến nghị để giúp nhà đầu tư hay các tổ chức có thể quản lý rủi ro hiệu quả hay nhìn thấy cơ hội đầu tư khi các thị trường có biến động mạnh.

1.2.2 Mục tiêu chi tiết

• Đánh giá cấu trúc phụ thuộc giữa TTCK Việt Nam và TTCK Nhật Bản nhằm xác định mức độ phản ứng của TTCK Việt Nam khi TTCK Nhật Bản biến động mạnh.

• So sánh sự khác nhau và kiểm định mức độ tác động của TTCK Nhật Bản đến TTCK Việt Nam ở các giai đoạn trước, trong và sau đại dịch Covid-19 và cuộc xung đột Nga - Ukraina.

• Đưa ra các kiến nghị nhằm giúp nhà đầu tư hay tổ chức quản lý các danh mục, rủi ro và cơ hội đầu tư.

1.3 Đối tượng nghiên cứu

Đối tượng nghiên cứu của bài tiểu luận là chuỗi tỷ suất lợi nhuận của chỉ số thị trường chứng khoán Việt Nam (VNINDEX) và thị trường chứng khoán Nhật Bản (Nikkei)

1.4 Phạm vi nghiên cứu

• Về không gian: Bài tiểu luận nghiên cứu mức độ phụ thuộc của TTCK Việt Nam khi TTCK Nhật Bản biến động.

• Về thời gian: Phạm vi thời gian được sử dụng để phân tích cấu trúc phụ thuộc là trong giai đoạn 04/2018 - 7/2024, dữ liệu được thu thập là dữ liệu theo ngày nhằm đánh giá khách quan hơn.

1.5 Kết cấu của đề tài

Chương 2: Cơ sở lý thuyết và tổng quan nghiên cứu

2.1 Một số lý thuyết liên quan

2.1.1 Lý thuyết về thị trường chứng khoán

Cũng như nhiều thị trường khác Thị trường chứng khoán cũng bao gồm một tập hợp bao gồm những người mua và người bán cổ phiếu (hay chứng khoán), thứ đại diện cho quyền sở hữu của họ đối với một doanh nghiệp; chúng có thể bao gồm các cổ phiếu được niêm yết trên sàn giao dịch chứng khoán đại chúng, hoặc những cổ phiếu được giao dịch một cách không công khai, ví dụ như cổ phần của một công ty tư nhân được bán cho các nhà đầu tư thông qua các nền tảng gọi vốn cộng đồng. Những khoản đầu tư trên thị trường chứng khoán hầu hết được thực hiện thông qua môi giới chứng khoán và nền tảng giao dịch điện tử.

Mục đích của thị trường chứng khoán là cung cấp một nơi mà bất kỳ ai cũng có thể mua và bán quyền sở hữu một phần trong một công ty đại chúng. Nó phân phối quyền kiểm soát một số công ty lớn cho hàng trăm triệu nhà đầu tư cá nhân. Và các quyết định mua và bán của những nhà đầu tư đó quyết định giá trị của các công ty đó. Thị trường chứng khoán có thể cho phép người mua và người bán thương lượng giá cả. Quá trình thương lượng này tối đa hóa sự công bằng cho cả hai bên bằng cách cung cấp cả giá bán cao nhất có thể và giá mua thấp nhất có thể tại một thời điểm nhất định. Mỗi sàn giao dịch theo dõi cung và cầu của các cổ phiếu được niêm yết tại đó.

Trên thực tế, thuật ngữ “thị trường chứng khoán” thường ám chỉ một trong những chỉ số thị trường chứng khoán chính, chẳng hạn như Chỉ số trung bình công nghiệp Dow Jones hoặc S&P 500 hay là đối với thị trường chứng khoán Việt Nam sẽ có chỉ số đại diện là VNINDEX. Những chỉ số này đại diện cho phần lớn thị trường chứng khoán và là thước đo sức khỏe của thị trường. Bởi vì khó có thể theo dõi từng công ty riêng lẻ, nên hiệu suất của các chỉ số này được coi là đại diện cho toàn bộ thị trường. Ta có thể thấy một tiêu đề tin tức nói rằng thị trường chứng khoán đã giảm hoặc thị trường chứng khoán đã đóng cửa tăng hoặc giảm trong ngày, điều này thường có nghĩa là các chỉ số thị trường chứng khoán đó đã tăng hoặc giảm và cổ phiếu trong chỉ số đã tăng hoặc giảm giá trị. Và các nhà đầu tư tham gia mua và bán cổ phiếu hy vọng sẽ kiếm được lợi nhuận thông qua biến động giá cổ phiếu này.

Dựa trên thời hạn luân chuyển vốn, thị trường chứng khoán được chia làm 2 loại đó là thị trường sơ cấp và thị trường thứ cấp.Thị trường sơ cấp là nơi chứng khoán được tạo ra. Tại thị trường này, các công ty bán (phát hành) cổ phiếu và trái phiếu mới cho công chúng lần đầu tiên, IPO (Initial Public Offering) là một ví dụ về thị trường sơ cấp. Các giao dịch này tạo cơ hội cho các nhà đầu tư mua chứng khoán từ ngân hàng đã thực hiện bảo lãnh phát hành ban đầu cho một cổ phiếu cụ thể. IPO xảy ra khi một công ty tư nhân phát hành cổ phiếu ra công chúng lần đầu tiên. Sau khi chứng khoán được phát hành và niêm yết cổ phiếu được tạo, các cổ phiếu và trái phiếu mới có thể được giao dịch trên thị trường thứ cấp thông qua các nhà môi giới hay nền tảng giao dịch điện tử.

2.1.2 Tổng quan thị trường chứng khoán (TTCK) Việt Nam

Tính đến nay, TTCK Việt Nam đã được hơn 25 năm thành lập và phát triển, dù đã trải qua nhiều thăng trầm song mức độ quan trọng của TTCK cho đến thời điểm hiện tại vẫn chưa có dấu hiệu suy giảm. Về thành quả, với thời gian dài hoạt động, TTCK Việt Nam đã khẳng định được vị thế của mình bằng những kết quả hết sức ấn tượng. Đầu tiên về quy mô, quy mô thị trường được mở rộng qua từng ngày, từ những năm đầu hình thành TTCK Việt Nam đã có những con số vô cùng khiêm tốn với chỉ 2 mã cổ phiếu được giao dịch lần lượt là của Công ty Cổ phần Cơ điện lạnh – REE và Công ty Cổ phần SAM Holdings – SAM, ở thời điểm đó mỗi tuần chỉ diễn ra 2 phiên giao dịch. Tuy nhiên, với sự ra đời của Trung tâm Giao dịch Chứng khoán Hà Nội (tiền thân của Sở Giao dịch Chứng khoán Hà Nội – HSX) vào ngày 8/3/2005, TTCK Việt Nam đã có những bước chuyển mình mạnh mẽ, cụ thể vào tháng 1/2006 cổ phiếu của Vinamilk lên sàn khiến giá trị vốn hóa của của HOSE tăng gấp đôi trong ngày.

Cũng trong năm này, 74 doanh nghiệp mới niêm yết trên sàn HOSE. Điều này giúp cho giá trị vốn hóa thị trường tăng từ 7.400 tỷ đồng lên 148.000 tỷ đồng. Chỉ số VN-Index cũng lên mức 752 điểm, tăng 144% chỉ sau một năm. Đỉnh điểm là vào ngày 12/3/2007, chỉ số VN-Index đạt đỉnh lịch sử khi kết thúc phiên với mức 1.170,69, con số này cao hơn khoảng 3.9 lần so với đầu năm 2006. Thế nhưng, 2 năm sau đó, cùng với sự khủng hoảng kinh tế toàn cầu, TTCK Việt Nam cũng theo đó mà bị ảnh hưởng, chỉ số VN-Index vào cuối năm 2008 đã giảm 66% so với cùng kỳ năm ngoái.

Sau giai đoạn khủng hoảng tài chính, thị trường bắt đầu nhịp tăng trưởng trở lại và chứng kiến nhiều sự thay đổi mới mẻ, nổi bật nhất là vào ngày 6/2/2012, chỉ số VN30 lần đầu tiên được ra mắt. VN30 là 30 cổ phiếu có giá trị vốn hóa lớn được niêm yết trên sàn HOSE, điều này đã đánh dấu bước chuyển mình theo tiêu chuẩn quốc tế, đồng thời giúp thị trường tăng thêm tính thanh khoản và minh bạch. Tiếp đến năm 2017, một hình thức chứng khoán mới là chứng khoán phái sinh được ra đời. Một loạt thay đổi trong thời gian trước cũng tạo ra những thay đổi tích cực. Cũng trong năm này, lịch sử chứng khoán Việt Nam chứng kiến hàng loạt con số liên tục “lập đỉnh”, VN-Index tăng 48% lên 984 điểm. Vốn hóa thị trường tăng hơn 70%, tỷ lệ trên GDP lần đầu vượt 50%. Tổng giá trị danh mục nhà đầu tư nước ngoài hơn 31 tỷ USD. Tiếp tục đỉnh cao, sau 10 năm, VN-Index lại thiết lập đỉnh mới vào ngày 9/4/2015 với 1.204 điểm. Đến năm 2023, tổng vốn hóa TTCK Việt Nam đến hết năm đạt hơn 240 tỷ USD, tương đương 56,4% GDP, chỉ số VN-Index chốt ở 1.129,93 điểm, cùng với đó là có khoảng 1900 công ty được niêm yết trên TTCK Việt Nam.

Thành quả kế tiếp có thể kể đến đó là cấu trúc của thị trường đang liên tục được hoàn thiện, việc áp dụng công nghệ số vào hệ thống giao dịch giúp hiệu suất và tính an toàn của các giao dịch được tăng lên đáng kể, cùng với đó là các quy định công bố thông tin chặt chẽ giúp cải thiện tính minh bạch đồng thời tăng thêm niềm tin cho các nhà đầu tư, thêm nữa TTCK Việt Nam đã xuất hiện thêm những loại sản phẩm mới ngoài sản phẩm cốt lõi là cổ phiếu như chứng khoán phái sinh, ETF và trái phiếu . Những thành quả này đã và đang giúp TTCK Việt Nam trở thành kênh huy động vốn trung và dài hạn vô cùng quan trọng cho các doanh nghiệp cũng như nền kinh tế Việt Nam.

2.1.3 Tổng quan thị trường chứng khoán (TTCK) Nhật Bản

Thị trường chứng khoán Tokyo (TSE) là một trong những sàn giao dịch chứng khoán lớn nhất và lâu đời nhất châu Á, đóng vai trò trung tâm trong nền kinh tế Nhật Bản. Được thành lập vào năm 1878, TSE đã chứng kiến nhiều giai đoạn thăng trầm cùng với sự phát triển của nền kinh tế Nhật Bản. TSE là nơi niêm yết của các tập đoàn lớn của Nhật Bản như Toyota, Sony, Mitsubishi, và nhiều công ty đa quốc gia khác.

TSE cung cấp một loạt các sản phẩm tài chính, bao gồm cổ phiếu, trái phiếu, quỹ đầu tư, và các sản phẩm phái sinh. Sàn giao dịch này cũng nổi tiếng với các chỉ số chứng khoán như Nikkei 225, một trong những chỉ số chứng khoán được theo dõi rộng rãi nhất trên thế giới. Nikkei 225 phản ánh giá trị thị trường của 225 công ty niêm yết lớn nhất trên TSE, và thường được xem là thước đo sức khỏe của nền kinh tế Nhật Bản.

> > Giống như nhiều thị trường chứng khoán khác, TSE cũng đã trải qua nhiều cuộc khủng hoảng và biến động. Tuy nhiên, với lịch sử lâu đời và cơ sở hạ tầng hiện đại, TSE vẫn giữ được vị thế là một trong những trung tâm tài chính quan trọng nhất châu Á. Trong những năm gần đây, TSE đã tích cực cải cách và hiện đại hóa các hệ thống giao dịch, nhằm tăng cường tính thanh khoản và hấp dẫn các nhà đầu tư nước ngoài >

2.1.4 Lý thuyết về cấu trúc phụ thuộc

Cấu trúc phụ thuộc đề cập đến cách mà các biến ngẫu nhiên hoặc điểm dữ liệu có liên quan hoặc phụ thuộc lẫn nhau. Cấu trúc phụ thuộc là một khái niệm chính trong thống kê, lý thuyết xác suất và nhiều lĩnh vực khác liên quan đến phân tích dữ liệu. Hiểu được cấu trúc phụ thuộc giữa các biến là tiền đề quan trọng để lập mô hình, phân tích và diễn giải dữ liệu một cách chính xác. Hầu hết các cấu trúc tương quan thường tồn tại dưới hai dạng chủ yếu là tuyến tính và phi tuyến tính.

Cấu trúc phụ thuộc tuyến tính đề cập đến mối quan hệ giữa hai hoặc nhiều biến trong đó sự phụ thuộc của chúng có thể được mô tả đầy đủ bằng mô hình tuyến tính. Điều này có nghĩa là những thay đổi trong một biến có thể được giải thích bằng những thay đổi tỷ lệ trong biến khác. Thông thường, để đo lường mức độ phụ thuộc tuyến tính thì hệ số tương quan Pearson được sử dụng rộng rãi. Pearson với giả định các chuỗi lợi suất của các tài sản tài chính có phân phối chuẩn hay phân phối Gauss và mối quan hệ giữa các tài sản tài chính là tuyến tính.

Phương pháp này đo lường tương quan song phương trong thời kỳ ổn định và kiểm định so sánh với tương quan trong và sau thời kỳ khủng hoảng. Nếu hệ số tương quan giai đoạn nghiên cứu sau tăng đáng kể so với giai đoạn trước đó, chứng tỏ tồn tại sự lây lan hay tồn tại sự đồng biến. Tuy nhiên, như đã nói ở trên các tỷ suất lợi nhuận trên thực tế thường không có phân phối chuẩn nên giả định của hệ số tương quan Pearson sẽ không phù hợp để giải thích cho cấu trúc phụ thuộc của các chỉ số này.

Cấu trúc phụ thuộc phi tuyến tính mô tả mối quan hệ giữa các biến mà mô hình tuyến tính không thể nắm bắt đầy đủ. Những mối quan hệ này có thể bao gồm các tương tác phức tạp hơn, có thể được xác định và mô hình hóa bằng nhiều kỹ thuật và công cụ khác nhau. Trong phạm vi bài tiểu luận này, chúng ta sẽ chỉ bàn luận về mô hình Copula. Copula là các hàm nắm bắt cấu trúc phụ thuộc giữa các biến ngẫu nhiên, cho phép mô hình hóa các mối phụ thuộc phức tạp vượt ra ngoài mối tương quan tuyến tính.

2.2 Lý thuyết về phương pháp Copula có điều kiện

2.2.1 Tổng quan về phương pháp Copula

Hàm copula đóng vai trò quan trọng trong lý thuyết xác suất và thống kê, giúp mô hình hóa quan hệ giữa các biến ngẫu nhiên đa chiều và xác định phân bố chung của chúng. Nó là công cụ thiết yếu trong phân tích đa biến, cho phép xác định mối quan hệ giữa các biến ngẫu nhiên một cách linh hoạt và hiệu quả. Copula dần được các nhà nghiên cứu phát triển theo đa dạng phương pháp khác nhau để linh hoạt với mục đích nghiên cứu, tiêu biểu là mô hình GARCH – copula, phát triển từ những năm 2000, dựa trên các mô hình ARCH và GARCH để mô hình hóa sự biến động của các biến số tài chính. Bollerslev (1986) đã đề xuất mô hình GARCH để dự báo độ biến động cho các chuỗi thời gian tài chính, được xây dựng từ mô hình của Engle (1982). Mô hình GARCH thường được sử dụng để dự báo ngắn hạn, đòi hỏi cập nhật thường xuyên. Để mô hình hóa mối liên kết giữa các chuỗi thời gian tài chính, nhiều mô hình GARCH đa biến đã được phát triển, chia thành hai hướng tiếp cận chính: phân tích ma trận hiệp phương sai trực tiếp và qua biến trung gian. Các mô hình GARCH mở rộng như ARMA-GARCH, APGARCH, EGARCH, TGARCH, FIGARCH thuộc hướng tiếp cận đầu tiên. Hướng thứ hai bao gồm các mô hình như DCC-GARCH, CCC-GARCH.

Sklar (1959) lần đầu tiên giới thiệu hàm copula, xác lập nền tảng lý thuyết copula, với tên gọi “copula” xuất phát từ tiếng Latinh có nghĩa là “liên kết” hoặc “kết nối”. Schweizer và Wol (1981) là những người đầu tiên liên kết copula với sự phụ thuộc của các biến ngẫu nhiên. Copula là một hàm phân phối tích lũy đa biến, xác định trên không gian n chiều với mọi phân phối biên đều nằm trong khoảng [0,1], dùng để mô hình hóa sự phụ thuộc giữa các biến ngẫu nhiên. Hàm copula cho phép xây dựng phân phối đa biến với các phân phối biên khác nhau và mô tả cấu trúc phụ thuộc giữa chúng. Sự phát triển lý thuyết copula đã được thúc đẩy bởi sự tiến bộ trong tính toán và phần mềm hỗ trợ. Copula đã trở thành công cụ phổ biến trong việc xác định sự phụ thuộc giữa các tài sản tài chính cũng như giữa tài sản tài chính và tài sản thực. Các copula đối xứng thuộc họ Elip, như Copula Gaussian và Copula Student-t, được sử dụng phổ biến trong thị trường tài chính.

Tóm lại, copula là một công cụ mạnh mẽ trong mô hình hóa sự phụ thuộc giữa các biến ngẫu nhiên. Khả năng của copula trong việc biểu diễn sự phụ thuộc thông qua các tham số như 𝜃trong Frank copula cho phép chúng ta phân tích và hiểu rõ hơn về mối quan hệ giữa các biến. Các công cụ tính toán như dCopula() trong R giúp các nhà thống kê dễ dàng triển khai và áp dụng copula vào các bài toán thực tế, từ đó mang lại những kết quả chính xác và có giá trị trong nghiên cứu và ứng dụng. Dưới đây là một số biểu đồ hàm mật độ xác suất của Frank Copula:

library(copula)
d <- 2 
theta <- -9
fc <- frankCopula(theta, dim = d)
wireframe2(fc, FUN = pCopula, draw.4.pCoplines = FALSE, main = "Biểu đồ Wireframe của Hàm Phân Phối Tích Lũy (CDF) của Frank Copula")

wireframe2(fc, FUN = dCopula, delta = 0.001, main = "Biểu đồ Wireframe của Hàm Mật Độ Xác Suất (PDF) của Frank Copula")

contourplot2(fc, FUN = pCopula, main = "Biểu đồ Contour của Hàm Phân Phối Tích Lũy (CDF) của Frank Copula")

contourplot2(fc, FUN = dCopula, n.grip = 72, lwd = 1/2,  main = "Biểu đồ Contour của Hàm Mật Độ Xác Suất (PDF) của Frank Copula")

>Bốn biểu đồ khác nhau, mỗi biểu đồ thể hiện một khía cạnh của hàm mật độ xác suất. Hai biểu đồ trên cùng là biểu đồ khung dây (wireframe) và biểu đồ đường đồng mức (contour plot) của hàm mật độ xác suất với d=2 và θ=−9. Biểu đồ khung dây (trên cùng bên trái) cho thấy một mô hình ba chiều của hàm mật độ xác suất, trong khi biểu đồ đường đồng mức (trên cùng bên phải) thể hiện các đường đồng mức của hàm mật độ này. Hai biểu đồ dưới cùng là các biểu đồ đường đồng mức của hàm mật độ xác suất với các giá trị khác nhau của θ. Biểu đồ dưới cùng bên trái thể hiện hàm mật độ xác suất với θ=0, trong khi biểu đồ dưới cùng bên phải thể hiện hàm mật độ xác suất với θ =−9. >

set.seed(1946)
n <- 1000
U <- rCopula(n, copula = fc)
U0 <- rCopula(n, copula = setTheta(fc, value = 0 ))
U9 <- rCopula(n, copula = setTheta(fc, value = 9))
plot(U, xlab = quote(U[1]), ylab = quote(U[2]),main = "Biểu đồ phân tán của mẫu U với θ = -9")

plot(U0, xlab = quote(U[1]), ylab = quote(U[2]),main = "Biểu đồ phân tán của mẫu U0 với θ = 0 (copula độc lập)")

plot(U9, xlab = quote(U[1]), ylab = quote(U[2]),main = "Biểu đồ phân tán của mẫu U9 với θ = 9")

2.2.2 Các họ hàm cơ bản của Copula

Copula được chia thành nhiều họ hàm khác nhau tương ứng cho mỗi mục đích và dữ liệu phù hợp, các họ hàm của Copula cơ bản được chia thành 2 nhánh chính là: Copula họ Elip và họ Archimedean.

Copula được chia thành nhiều họ hàm khác nhau tương ứng cho mỗi mục đích và dữ liệu phù hợp, các họ hàm của Copula cơ bản được chia thành 2 nhánh chính là: Copula họ Elip và họ Archimedean.

Copula họ Elip

Copula họ Elip gồm Copula Gaussian và Copula Student-t, là những công cụ quan trọng trong phân tích tài chính. Copula Gaussian dựa trên phân phối chuẩn và thể hiện mối quan hệ tuyến tính giữa các biến ngẫu nhiên, phù hợp với các điều kiện thị trường bình thường. Khi các biến ngẫu nhiên gần giá trị trung bình, copula này cho thấy chúng có thể hoạt động độc lập, phản ánh sự tương quan khi thị trường ổn định. Trong khi đó, Copula Student-t, dựa trên phân phối Student-t, cho phép mô hình hóa các sự kiện ở phần đuôi của phân phối. Copula này đặc biệt hữu ích khi thị trường biến động mạnh, vì nó cho thấy sự phụ thuộc lớn hơn khi thị trường thay đổi nhiều, nhất là khi số bậc tự do nhỏ. Điều này làm cho Copula Student-t trở thành công cụ quan trọng trong việc mô hình hóa rủi ro tài chính, khi các sự kiện cực đoan có thể xảy ra.

Ta có thể đánh giá ý nghĩa của Copula họ Elip thông qua các biểu đồ bên dưới dựa trên dữ liệu được mô phỏng từ hàm Copula tương ứng.

library(knitr)
library(copula)
library(ggplot2)
library(plotly)
library(gridExtra)
library(VC2copula)
library(VineCopula)
library(gridGraphics)
library(png)
#Tạo hàm vẽ biểu đồ
scatter_plot <- function(random_data, cl) {
  ggplot(data.frame(random_data), aes(random_data[,1], random_data[,2])) +
  geom_point(alpha = 0.5, col = cl) +
  theme_minimal() +
  labs(x = "u", y = "v")
  } #for scatter

persp_plot <- function(copula_obj, file_name, cl) {
  png(file_name)
  persp(copula_obj, dCopula, 
        xlab = 'u', ylab = 'v', col = cl, ltheta = 120,  
        ticktype = "detailed", cex.axis = 0.8)
  dev.off()
  rasterGrob(readPNG(file_name),interpolate = TRUE)
}# for pdf

set.seed(123)
#Mô phỏng copula gauss với p=0.8
cop_nor <- normalCopula(param = 0.8, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_nor <- rCopula(copula = cop_nor,n = 7600)
#Mô phỏng copula với p=0.8
cop_std <- tCopula(param = 0.8, dim = 2, df = 1)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_std <- rCopula(copula = cop_std,n = 7600)

#Vẽ biểu đồ 
png('nors.png') #Lưu ảnh biểu đồ
scatter_plot(random_nor, '#99FFCC')
dev.off()
nors <- rasterGrob(readPNG('nors.png'), interpolate = TRUE) #Chuyển ảnh sang dạng Grob

png('stds.png')
scatter_plot(random_std, '#CCFF33')
dev.off()
stds <- rasterGrob(readPNG('stds.png'), interpolate = TRUE)

nor_per <- persp_plot(cop_nor, 'norp.png', '#99FFCC')
std_per <- persp_plot(cop_std, 'stdp.png', '#CCFF33')

legend <- legendGrob(
  labels = c("Gauss", "Student"), pch = 15,
  gp = gpar(col = c('#99FFCC', '#CCFF33'), fill = c('#99FFCC', '#CCFF33'))
)

grid.arrange(nors, nor_per, stds, std_per, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 1: Biểu đồ phân tán và phối cảnh PDF của Copula họ Elip", 
                 gp = gpar(fontsize = 15, font = 2))
             )

Copula họ Archimedean

Copula họ Archimedean bao gồm nhiều hàm như Copula Clayton, Gumbel, và Joe, mỗi loại thể hiện một kiểu phụ thuộc đặc biệt giữa các biến ngẫu nhiên. Copula Clayton tập trung vào sự phụ thuộc mạnh mẽ ở phần đuôi dưới của phân phối, làm cho nó phù hợp khi mô hình hóa mối quan hệ giảm giá mạnh giữa các biến, hữu ích trong các tình huống thị trường giảm. Ngược lại, Copula Gumbel thể hiện sự phụ thuộc mạnh mẽ ở phần đuôi trên, thích hợp khi một biến tăng giá có thể kéo theo biến khác cũng tăng giá mạnh, thường được sử dụng trong thị trường có xu hướng tăng. Copula Joe, tương tự như Gumbel, cũng tập trung vào phần đuôi trên, mô tả mối quan hệ mạnh mẽ khi các biến tăng giá đồng thời. Các copula khác như Frank, Rot-Clayton, và các copula phức tạp như BB1, BB6, BB7, BB8, kết hợp các đặc điểm của nhiều copula cơ bản, cho phép mô hình hóa các mối quan hệ phức tạp và đa dạng hơn giữa các biến ngẫu nhiên trong tài chính, giúp hiểu rõ hơn về sự tương tác và rủi ro trong môi trường tài chính biến đổi.

Để có thể đánh giá ý nghĩa của từng họ Copula riêng biệt một cách trực quan, em tiến hành vẽ biểu đồ dựa trên dữ liệu được mô phỏng của mỗi họ Copula tương ứng.

set.seed(123)
#Mô phỏng copula clayton với p=4
cop_clay <- claytonCopula(param = 4, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_clay <- rCopula(copula = cop_clay,n = 7600)
#Mô phỏng copula gumbel với p=5
cop_gum <- gumbelCopula(param = 5, dim = 2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_gum <- rCopula(copula = cop_gum,n = 7600)

#Vẽ biểu đồ 

png('clays.png')
scatter_plot(random_clay,'#CC33FF')
dev.off()
clays <- rasterGrob(readPNG('clays.png'), interpolate = TRUE)

png('gums.png')
scatter_plot(random_gum,'#FF6699')
dev.off()
gums <- rasterGrob(readPNG('gums.png'), interpolate = TRUE)

clayp <- persp_plot(cop_clay,'clayp.png','#CC33FF')

gump <- persp_plot(cop_gum,'gump.png','#FF6699')


legend <- legendGrob(
  labels = c("Clayton", "Gumbel"), pch = 15,
  gp = gpar(col = c('#CC33FF', '#FF6699'), fill = c('#CC33FF', '#FF6699'))
)

grid.arrange(clays, clayp, gums, gump, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3, 4, 5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 2: Biểu đồ phân tán và PDF của Copula Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula survival gumbel
cop_surgum <- VC2copula::surGumbelCopula(param = 5)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_surgum <- rCopula(copula = cop_surgum,n = 7600)
#Mô phỏng copula survival clayton  với p=4
cop_surclay <- VC2copula::surClaytonCopula(param = 4)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_surclay <- rCopula(copula = cop_surclay,n = 7600)

png('surclays.png')
scatter_plot(random_surclay,'#FFE4B5')
dev.off()
surclays <- rasterGrob(readPNG('surclays.png'), interpolate = TRUE)

png('surgums.png')
scatter_plot(random_surgum,'#0000EE')
dev.off()
surgums <- rasterGrob(readPNG('surgums.png'), interpolate = TRUE)

surclayp <- persp_plot(cop_surclay, "surclayp.png","#FFE4B5")
surgump <- persp_plot(cop_surgum, "surgump.png","#0000EE")

legend <- legendGrob(labels = c("Survival Clayton","Survival Gumbel"), pch = 15,
                    gp = gpar(col = c('#FFE4B5','#0000EE'), fill = c('#FFE4B5','#0000EE' )))

grid.arrange(surclays, surclayp,surgums, surgump, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 3: Biểu đồ phân tán và PDF của Copula Sur Clayton và Gumbel", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula Frank
cop_frank <- frankCopula(param = 9.2)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_frank <- rCopula(copula = cop_frank,n = 7600)
#Mô phỏng copula survival clayton  với p=4
cop_joe <- joeCopula(param = 3)
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_joe <- rCopula(copula = cop_joe,n = 7600)

png('franks.png')
scatter_plot(random_frank,'#33CC33')
dev.off()
franks <- rasterGrob(readPNG('franks.png'), interpolate = TRUE)

png('joes.png')
scatter_plot(random_joe,'#FFCC66')
dev.off()
joes <- rasterGrob(readPNG('joes.png'), interpolate = TRUE)

frankp <- persp_plot(cop_frank, "frankp.png","#33CC33")
joep <- persp_plot(cop_joe, "joep.png","#FFCC66")

legend <- legendGrob(labels = c("Franl","Joe"), pch = 15,
                    gp = gpar(col = c('#33CC33','#FFCC66'), fill = c('#33CC33','#FFCC66')))

grid.arrange(franks, frankp,joes, joep, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 4: Biểu đồ phân tán và  PDF của Copula Frank và Joe", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula BB1
cop_bb1 <- VC2copula::BB1Copula(param = c(2,1.7))
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_bb1 <- rCopula(copula = cop_bb1,n = 7600)
#Mô phỏng copula BB6
cop_bb6 <- VC2copula::BB6Copula(param = c(2,4.5))
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_bb6 <- rCopula(copula = cop_bb6,n = 7600)

png('bb1s.png')
scatter_plot(random_bb1,'#FF6600')
dev.off()
bb1s <- rasterGrob(readPNG('bb1s.png'), interpolate = TRUE)

png('bb6s.png')
scatter_plot(random_bb6,'#999900')
dev.off()
bb6s <- rasterGrob(readPNG('bb6s.png'), interpolate = TRUE)

bb1p <- persp_plot(cop_bb1, "bb1p.png","#FF6600")
bb6p <- persp_plot(cop_bb6, "bb6p.png","#999900")

legend <- legendGrob(labels = c("BB1","BB6"), pch = 15,
                    gp = gpar(col = c('#FF6600','#999900'), fill = c('#FF6600','#999900')))

grid.arrange(bb1s, bb1p,bb6s, bb6p, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 5: Biểu đồ phân tán và  PDF của Copula BB1 và BB6", 
                 gp = gpar(fontsize = 15, font = 2))
             )

set.seed(123)
#Mô phỏng copula BB7
cop_bb7 <- VC2copula::BB7Copula(param = c(2,4.5))
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_bb7 <- rCopula(copula = cop_bb7,n = 7600)
#Mô phỏng copula BB8
cop_bb8 <- VC2copula::BB8Copula(param = c(4,0.8))
#Mô phỏng 7600 quan sát ngẫu nhiên dựa trên copula có sẵn
random_bb8 <- rCopula(copula = cop_bb8,n = 7600)

png('bb7s.png')
scatter_plot(random_bb7,'#FF0000')
dev.off()
bb7s <- rasterGrob(readPNG('bb7s.png'), interpolate = TRUE)

png('bb8s.png')
scatter_plot(random_bb8,'#33CCCC')
dev.off()
bb8s <- rasterGrob(readPNG('bb8s.png'), interpolate = TRUE)

bb7p <- persp_plot(cop_bb7, "bb7p.png","#FF0000")
bb8p <- persp_plot(cop_bb8, "bb8p.png","#33CCCC")

legend <- legendGrob(labels = c("BB7","BB8"), pch = 15,
                    gp = gpar(col = c('#FF0000','#33CCCC'), fill = c('#FF0000','#33CCCC')))

grid.arrange(bb7s, bb7p,bb8s, bb8p, legend, ncol = 3,
             layout_matrix = rbind(c(1,2,5),c(3,4,5)),
             widths = c(2,2,1),
             top = textGrob("Hình 6: Biểu đồ phân tán và  PDF của Copula BB7 và BB8", 
                 gp = gpar(fontsize = 15, font = 2))
             )

d <- 2 
theta <- -9
fc <- frankCopula(theta, dim = d)
wireframe2(fc, FUN = pCopula, draw.4.pCoplines = FALSE, main = "Biểu đồ Wireframe của Hàm Phân Phối Tích Lũy (CDF) của Frank Copula")

wireframe2(fc, FUN = dCopula, delta = 0.001, main = "Biểu đồ Wireframe của Hàm Mật Độ Xác Suất (PDF) của Frank Copula")

contourplot2(fc, FUN = pCopula, main = "Biểu đồ Contour của Hàm Phân Phối Tích Lũy (CDF) của Frank Copula")

contourplot2(fc, FUN = dCopula, n.grip = 72, lwd = 1/2,  main = "Biểu đồ Contour của Hàm Mật Độ Xác Suất (PDF) của Frank Copula")

2.2.3 Phụ thuộc đuôi

Cấu trúc phụ thuộc đuôi (tail dependence) mô tả mức độ mà các biến ngẫu nhiên có xu hướng cùng đạt các giá trị cực đoan cụ thể là ở đuôi góc phần tư phía trên bên phải hoặc ở phần đuôi góc phần tư phía dưới bên trái của phân phối hai biến, chẳng hạn như cùng tăng mạnh hoặc giảm mạnh, tại các phần đuôi của phân phối xác suất của chúng. Cấu trúc phụ thuộc đuôi giúp các nhà phân tích và nhà đầu tư hiểu rõ hơn về rủi ro cực đoan, như các sự kiện tài chính bất thường có thể xảy ra đồng thời, chẳng hạn như sự sụt giảm đồng loạt của các tài sản trong danh mục đầu tư. Có hai loại phụ thuộc đuôi chính:

• Phụ thuộc đuôi trên (Upper tail dependence): Đo lường mức độ các biến ngẫu nhiên có xu hướng đạt các giá trị cực đại đồng thời. Điều này quan trọng khi xem xét các kịch bản mà các tài sản tài chính cùng tăng giá mạnh, có thể dẫn đến sự tích lũy rủi ro khi thị trường tăng mạnh.

• Phụ thuộc đuôi dưới (Lower tail dependence): Đo lường mức độ các biến ngẫu nhiên có xu hướng đạt các giá trị cực tiểu đồng thời. Điều này đặc biệt quan trọng trong việc quản lý rủi ro tài chính, vì sự sụt giảm giá trị đồng loạt của các tài sản có thể gây ra thua lỗ lớn cho nhà đầu tư.

Các hệ số phụ thuộc đuôi, thường ký hiệu là \(\lambda_U\) cho phụ thuộc đuôi trên và \(\lambda_L\) cho phụ thuộc đuôi dưới, được định nghĩa như sau:

• Hệ số phụ thuộc đuôi trên: \[ \lambda_U = \lim_{u \to 1} P(Z_1 > F_1^{-1}(u) \mid Z_2 > F_2^{-1}(u)) \]

• Hệ số phụ thuộc đuôi dưới: \[ \lambda_L = \lim_{u \to 1} P(Z_1 < F_1^{-1}(u) \mid Z_2 < F_2^{-1}(u)) \]

Trong đó, \(Z_1\) và \(Z_2\) là các biến ngẫu nhiên liên tục, \(F_1\) và \(F_2\) là các hàm phân phối biên tương ứng. Hệ số \(\lambda_U\) hoặc \(\lambda_L\) bằng 0 cho thấy không có phụ thuộc đuôi, trong khi giá trị trong khoảng (0,1] cho thấy mức độ phụ thuộc đuôi tồn tại. Sự phụ thuộc đuôi là một thuộc tính quan trọng của copula, vì nó là bất biến dưới các phép biến đổi tăng nghiêm ngặt của các biến ngẫu nhiên. Điều này có nghĩa là, ngay cả khi các biến ngẫu nhiên được biến đổi theo cách không làm thay đổi thứ tự của chúng, sự phụ thuộc đuôi vẫn không thay đổi.

2.3 Tổng quan nghiên cứu nghiên cứu liên quan

2.3.1 Tổng quan nghiên cứu nước ngoài

Dự báo là nghệ thuật tiên đoán nhưng mang nhiều yếu tố khoa học, nghiên cứu trên cơ sở xử lý logic các dữ liệu, số liệu được thu thập liên tục trong quá khứ, hiện tại và tương lai để dự đoán các xu hướng vận động của các đối tượng trong tương lai nhờ vào một số mô hình toán học (Định lượng). Tuy nhiên, dự báo cũng có thể là một dự đoán chủ quan, có trực giá về tương lai (Định tính) của người dự báo, người ta thường cố loại trừ nó đi.Không có cách nào để xác định tương lai là một cách chắc chắn, nhưng dưới dòng chảy lịch sử, con người ta đã nghiên cứu, khắc phục và nghiên cứu liên tục nhiều mô hình để khắc phục các “khuyết tật” trước đó của bản thân nó – Mô hình chuỗi thời gian.

Mô hình chuỗi thời gian thực sự rất hữu ích trong các lĩnh vực của đời sống Tuy nhiên, để khắc phục những hạn chế của các mô hình này, nhiều nhà nghiên cứu đã tiến hành các cải tiến và mở rộng nhằm nâng cao độ chính xác và hiệu quả của dự báo Nghiên cứu về biến động điều kiện tự hồi quy (ARCH) và các mô hình liên quan có thể kể đến như (ARIMA) đã nhận được sự quan tâm rộng rãi từ cộng đồng khoa học tài chính quốc tế. Những đóng góp của các nhà nghiên cứu trong lĩnh vực này đã mở ra nhiều hướng đi mới và cải thiện đáng kể hiểu biết của chúng ta về các hiện tượng tài chính phức tạp.

Tim Bollerslev là một trong những nhà nghiên cứu tiên phong trong lĩnh vực này. Các công trình của ông, như “A Conditionally Heteroskedastic Time Series Model for Speculative Prices and Rates of Return” (Bollerslev, Aconditionally Heteroskedastic Time Series Model for Speculative, 1990) và (Bollerslev, Modelling The Coherence in Short-run Nominal Exchange Rates A Multivariate Generalized ARCH Model, 1990), đã giới thiệu các mô hình chuỗi thời gian tiên tiến để nắm bắt sự phụ thuộc trong dữ liệu tài chính. Bollerslev cũng đã đóng góp quan trọng vào việc hiểu rõ cấu trúc tương quan trong các mô hình biến động điều kiện tự hồi quy tổng quát (1988). Wai Keung Li cũng đã có những đóng góp đáng kể trong việc nghiên cứu các mô hình chuỗi thời gian phi tuyến với biến động điều kiện (Li, 1994). Các bài báo của ông như “On Fractionally Integrated Autoregressive Moving-Average Time Series Models with Conditional Heteroscedasticity” (1997), đã khám phá các tính chất và phân phối của các tự tương quan dư bình phương trong các mô hình này. Công trình của Li cung cấp các kỹ thuật ước lượng quan trọng cho các quá trình đơn vị gốc với sai số GARCH. Kế thừa những nghiên cứu đó, Mika Meitz, cùng với Pentti Saikkonen, đã tập trung vào ước lượng tham số và phân tích ổn định trong các mô hình tự hồi quy phi tuyến với sai số biến động điều kiện (Mika Meitz,2008). Các bài báo “Parameter Estimation in Nonlinear AR-GARCH Models” (2008) và “Stability of Nonlinear AR‐GARCH Models” (2008) của họ đã đóng góp vào lý thuyết ước lượng tiệm cận cho các mô hình này, cung cấp những hiểu biết sâu sắc về tính phi tuyến và biến động điều kiện trong dữ liệu chuỗi thời gian.

Sở dĩ chúng ta gọi các mô hình này là chuỗi thời gian (time series) là bởi mô hình được áp dụng các yếu tố đặc thù có liên quan đến thời gian Một mô hình chuỗi thời gian thường dự báo dựa trên một giả định rằng các quy luật trong quá khứ sẽ được lặp lại trong tương lai . Do đó xây dựng mô hình chuỗi thời gian là chúng ta đang mô hình hóa mối quan hệ trong quá khứ giữa biến độc lập (biến đầu vào) và biến phụ thuộc (biến mục tiêu). Dựa vào mối quan hệ này để dự đoán giá trị trong tương lai của biến phụ thuộc. Mô hình ở đây tôi đề xuất, có thể đơn giản trong cách vận hành, sử dụng cũng như mạnh mẽ trong việc dụ báo các giá trị trong tương lai – ARIMA, trước đó đã có nhiều nghiên cứu, ứng dụng của mô hình vào các khía cạnh khác nhau của đời sống. Cụ thể, Mehdi Khashei đã cùng với các đồng tác giả tiến hành nghiên cứu tích hợp các mô hình ARIMA với Mạng Nơ-ron Nhân tạo (ANN) và logic Mờ. Các công trình của họ, bao gồm các nghiên cứu được công bố vào các năm 2009, 2013 và 2015, tập trung vào việc khắc phục các hạn chế tuyến tính và dữ liệu của mô hình ARIMA. Bằng cách kết hợp ANN và logic Mờ, họ đã cải thiện đáng kể độ chính xác của các dự báo trong thị trường tài chính. Không chỉ với lĩnh vực tài chính, ứng dụng mô hình trong dự đoán các biến động thực tế cũng mang lại những kết quả khả quan. Mohammad Ebrahim Banihabib và các đồng nghiệp của ông đã đề xuất các mô hình tự hồi quy tuyến tính và phi tuyến mở rộng để dự báo tiêu thụ nước đô thị và lưu lượng đỉnh hàng ngày. Nghiên cứu của họ sử dụng các phương pháp ARIMA và Nonlinear Auto-Regressive Exogenous (NARX) để nâng cao độ chính xác của dự báo. Các công trình của họ, được công bố vào năm 2019 và 2020, đã đóng góp quan trọng vào việc cải thiện dự báo trong các ứng dụng thực tế liên quan đến quản lý nước. Những nghiên cứu trên đây thực sự là những bước đột phá trong dự báo những “nước đi” của số liệu, giúp con người có thêm dữ kiện trong việc phân tích, đánh giá khách quan hơn một hiện tượng.

2.3.2. Tổng quan nghiên cứu trong nước

Nghiên cứu về mô hình dự báo chuỗi thời gian đã được áp dụng rộng rãi trong nhiều lĩnh vực tại Việt Nam, đặc biệt là trong dự báo kinh tế và nông nghiệp. Các nghiên cứu này chủ yếu sử dụng các mô hình như ARIMA (Auto-Regressive Integrated Moving Average), GARCH (Generalized AutoRegressive Conditional Heteroskedasticity), và các mô hình tích hợp với các phương pháp hiện đại như mạng nơron nhân tạo (ANN) và logic mờ (Fuzzy Logic). Tại Việt Nam, trong nhưng năm gần đây các công ty điện lực rất quan tâm đến dự báo điện năng. Tuy nhiên, phần lớn các công ty dựa vào kinh nghiệm tích lũy cũng như những phương pháp thống kê đơn giản như ngoại suy, hồi quy và hệ số đàn hội. Phương pháp mạng nơ ron nhân tạo cũng được bắt đầu sử dụng nhưng mới ở dạng mô hình đơn lẻ. Tuy nhiên, nếu nhắc đến nghiên cứu về mô hình dự báo, Việt Nam ta cũng có nhưng thành tự đáng kể. Phúc và cộng sự (2022) đã thực hiện một phân tích hồi quy meta để đánh giá sự biến đổi trong các ước lượng hiệu quả kỹ thuật (TE) trong các nghiên cứu về sản xuất lúa gạo (Phuc , 2022). Họ phân tích 443 quan sát từ 175 nghiên cứu được công bố trong ba thập kỷ. Nghiên cứu này xác định một số yếu tố ảnh hưởng đến ước lượng TE, bao gồm sự không đồng nhất của nghiên cứu, loại dữ liệu và các khu vực địa lý. Hien, Kawaguchi và Suzuki tập trung vào việc đánh giá hiệu quả kỹ thuật của sản xuất lúa gạo tại Đồng bằng sông Cửu Long, Việt Nam bằng phương pháp phân tích biên ngẫu nhiên (Hien, 2003). Nghiên cứu này nhằm xác định mức độ hiệu quả và các yếu tố ảnh hưởng đến hiệu quả trong số các nông dân trồng lúa ở khu vực này.

Về ứng dụng ARIMA để dự báo, tiêu biểu có Huỳnh Tấn Nguyên và Nguyễn Văn Lượng đã nghiên cứu và dự báo chỉ số giá tiêu dùng CPI trong quý 2 năm 2016. Kết quả thu về cho thấy việc ứng dụng mô hình hợp lý đã cho những tầm nhìn khả qua trong việc lập kế hoạch kinh doanh, kế hoạch dự trữ hàng tồn kho để phục vụ tốt nhất cho nhu cầu của người tiêu dùng, đối với các nhà đầu tư khi xác định được CPI có thể tính toán tỷ lệ sinh lời của dự án đầu tư một cách hợp lý. Phức tạp hơn, Hương và Trung (2015) đã giới thiệu mô hình ARIMA-GARCH hiện đang được sử dụng khá phổ biến trong việc dự báo các chuỗi thời gian bên cạnh các mô hình đơn lẻ là mô hình ARIMA và mô hình GARCH; đồng thời ứng dụng được mô hình này trong đúng kết, dự báo chỉ số VN-Index của thị trường chứng khoán Việt Nam trong thực tế. Kết quả trả về cho thấy mô hình kết hợp ARIMA(1,1,1,)-GARCH(1,1) cho kết quả dự báo tốt hơn hai mô hình đơn lẻ còn lại là ARIMA(1,1,1) và GARCH(1,1).

Chương 3: Thiết kế nghiên cứu

3.1 Dữ liệu nghiên cứu

Dữ liệu nghiên cứu được sử dụng trong bài tiểu luận là hai chuỗi tỷ suất lợi nhuận của TTCK Việt Nam được đại diện bởi chỉ số giá VNINDEX và TTCK Nhật Bản được đại diện bởi chỉ số giá NK225. Giai đoạn nghiên cứu là từ tháng 4 năm 2018 đến đầu tháng 7 năm 2024. Với mục đích đo lường mức độ tác động của các cú sốc kinh tế nên nghiên cứu được thiết kế chia thành 4 giai đoạn chính là giai đoạn trước các cú sốc kinh tế bắt đầu từ đầu 2018 đến cuối năm 2019, giai đoạn thứ hai là trong dịch Covid - 19 từ đầu 2020 đến cuối năm 2021, tiếp theo là giai đoạn xung đột Nga - Ukraine bắt đầu từ 2022 đến tháng 1/2023, cuối cùng là giai đoạn sau cú sốc tiếp nối từ tháng 1/2023 đến đầu tháng 7/2024. Giá của hai chỉ số chứng khoán được thu thập từ trang web (Investing.com). Tỷ suất lợi nhuận được tính theo công thức Rt = ln(Pt) - ln(Pt-1), với Pt là giá chỉ số chứng khoán ở thời điểm t.

3.2 Phương pháp nghiên cứu

Phương pháp nghiên cứu của bài luận được trình bày theo thứ tự với bốn bước cụ thể là: (1) Xác định mô hình phân phối biên, (2) Kiểm định tính phù hợp của mô hình phân phối biên, (3) Ước lượng tham số mô hình copula, (4) Lựa chọn mô hình copula phù hợp.

3.2.1 Xác định mô hình phân phối biên

Nghiên cứu của chúng ta đang thực hiện một loạt các bước để xác định mô hình biên ARMA – GJR – GARCH thích hợp nhất cho dữ liệu chuỗi lợi suất. Các bước cụ thể như sau:

• Kiểm tra sự tồn tại của hiệu ứng ARCH: Sử dụng kiểm định Lagrange Multiplier (LM) để xác định liệu dữ liệu có hiệu ứng biến động không đều hay không.

• Xác định bậc p và q cho mô hình ARMA: Dựa vào các tiêu chuẩn thông tin như AIC (Akaike Information Criterion) và BIC (Bayesian Information Criterion), chúng ta lựa chọn bậc p và q phù hợp để đặc trưng cho phần tự hồi quy (AR) và phần trung bình chuyển động (MA) của chuỗi lợi suất.

• Xây dựng mô hình ARMA (p,q) – GJR – GARCH (r,m): Với các bậc p và q đã xác định cho mô hình ARMA, mô hình này được mở rộng bằng cách thêm phần GJR – GARCH (r,m) để mô hình hóa hiệu ứng GARCH và hiệu ứng đòn bẩy (leverage effect).

• Xác định phân phối cho phần dư: Lựa chọn phân phối phù hợp cho phần dư của mô hình GJR – GARCH, bao gồm các phân phối như Normal, Student-t, Skewed Student-t, GED (Generalized Error Distribution), và sGED (Skewed Generalized Error Distribution).

• Lựa chọn mô hình biên phù hợp nhất: Sử dụng các tiêu chuẩn thông tin như AIC, BIC, SIC và HQIC để chọn ra mô hình biên ARMA – GJR – GARCH phù hợp nhất. Mô hình có giá trị tiêu chuẩn thông tin thấp nhất sẽ được coi là phù hợp nhất với dữ liệu.

• Mô hình biên ARMA – GJR – GARCH hoàn chỉnh: Sau khi xác định mô hình biên thích hợp nhất, chúng ta có một mô hình ARMA – GJR – GARCH hoàn chỉnh mô tả các đặc tính quan trọng của chuỗi lợi suất, bao gồm đuôi dày, đối xứng và bất đối xứng, cùng hiệu ứng đòn bẩy.

Quá trình này giúp chúng ta xác định mô hình biên phù hợp nhất để mô tả các đặc điểm quan trọng của dữ liệu chuỗi lợi suất.

3.2.2 Kiểm định tính phù hợp cho mô hình phân phối biên

Sau khi đã xác định mô hình biên tối ưu cho mỗi chuỗi lợi suất, chúng ta tiến hành tính toán phần dư chuẩn hóa từ mô hình này, ký hiệu là cặp dữ liệu \(z_{1t}\) và \(z_{2t}\) tại thời điểm t. Tiếp theo, sử dụng hàm phân phối biên thực nghiệm để chuyển đổi \(z_{1t}\) và \(z_{2t}\) thành giá trị xác suất hoặc giá trị tích phân: \(u_t = F_1(z_{1t} | \Omega_{t-1})\) và \(v_t = F_2(z_{2t} | \Omega_{t-1})\).

Với các giá trị xác suất \(u_t\) và \(v_t\), chúng ta xây dựng mô hình copula phi tham số \(C(u_t, v_t)\), mô hình này mô tả mức độ phụ thuộc giữa \(u_t\) và \(v_t\) mà không cần xác định các tham số cụ thể của copula. Mô hình copula phi tham số có tính linh hoạt cao hơn và phù hợp với nhiều tình huống so với mô hình copula tham số.

Gọi \(C\) là hàm copula, với \(u_t, v_t \in [0, 1]^2\) và các vi phân từng phần \(\frac{\partial C(u_t, v_t)}{\partial u_t}\), \(\frac{\partial C(u_t, v_t)}{\partial v_t}\) tồn tại cho từng biến \(u_t, v_t\): \[ 0 \leq \frac{\partial C(u_t, v_t)}{\partial u_t} \leq 1, \quad 0 \leq \frac{\partial C(u_t, v_t)}{\partial v_t} \leq 1 \]

Trong đó, \(\frac{\partial C(u_t, v_t)}{\partial u_t}\) và \(\frac{\partial C(u_t, v_t)}{\partial v_t}\) lần lượt là hàm phân phối biên của các biến \(u_t\) và \(v_t\).

Cuối cùng, để đảm bảo tính phù hợp của mô hình copula và hàm phân phối biên, chúng ta sử dụng các kiểm định thống kê bao gồm:

• Kiểm định Anderson – Darling (A-D): Kiểm tra sự khớp giữa phân phối biên dự kiến và phân phối biên thực nghiệm. Giá trị p lớn cho thấy không có bằng chứng đủ để bác bỏ tính phù hợp của hàm phân phối biên.

• Kiểm định Cramer – von Mises (Cv-M): Đo lường sự khác biệt giữa phân phối biên dự kiến và phân phối biên thực nghiệm. Giá trị p của kiểm định này cũng phản ánh mức độ phù hợp của hàm phân phối biên.

• Kiểm định Kolmogorov – Smirnov (K-S): Kiểm tra sự khác biệt giữa phân phối biên dự kiến và phân phối biên thực nghiệm. Tương tự như A-D và Cv-M, giá trị p của K-S cũng giúp đánh giá tính phù hợp của hàm phân phối biên.

3.2.3 Ước lượng các tham số

Ước lượng tham số cho mô hình copula có thể thực hiện bằng phương pháp suy luận cận biên (IFM), một kỹ thuật phổ biến cho phép phân tách việc ước lượng tham số của hàm phân phối biên và tham số của copula. Phương pháp này đặc biệt hữu ích cho các mô hình điều kiện phức tạp hơn. Ban đầu, chúng ta ước lượng các tham số của hàm phân phối biên cho từng biến độc lập, sử dụng các phương pháp như MLE hoặc phương pháp moments, tùy thuộc vào loại phân phối. Sau khi hoàn tất ước lượng các tham số của hàm phân phối biên, chúng ta xây dựng hàm suy luận cận biên (IFM).

Hàm IFM liên quan đến hàm mật độ xác suất (PDF) hoặc hàm khối lượng xác suất (PMF) của từng biến độc lập, dựa trên các tham số đã được ước lượng. Hàm IFM cho phép tính toán các đạo hàm của hàm copula theo các tham số của hàm phân phối biên. Khi có hàm IFM, chúng ta có thể áp dụng các phương pháp để ước lượng các tham số của copula, phân biệt từ tham số của hàm phân phối biên đã ước lượng trước đó. Phương pháp IFM hỗ trợ việc xử lý các mô hình điều kiện phức tạp bằng cách ước lượng riêng biệt cho các thành phần của copula và phân phối biên.

Quá trình IFM gồm hai bước chính: Bước 1: Ước lượng các tham số của hàm phân phối biên bằng phương pháp MLE:

• \[ \hat{\theta}_1 = \arg\max \sum_{t=1}^T \log f_1(z_{1t} \mid \Omega_{t-1}, \theta_1)\]

• \[ \hat{\theta}_2 = \arg\max \sum_{t=1}^T \log f_2(z_{2t} \mid \Omega_{t-1}, \theta_2)\]

Bước 2: Ước lượng tham số của copula

\[\hat{\theta}_c = \arg\max \sum_{t=1}^T (F_1(z_{1t} \mid \Omega_{t-1}; \hat{\theta}_1), F_2(z_{2t} \mid \Omega_{t-1}; \hat{\theta}_2), \theta_c)\]

Trong đó, \(l_c\) là hàm log-likelihood của hàm mật độ xác suất của copula \(c\), và \(l_1\), \(l_2\) là hàm log-likelihood của các phân phối biên \(F_1\) và \(F_2\).

3.2.4 Lựa chọn mô hình Copula phù hợp

Để lựa chọn mô hình copula phù hợp, việc sử dụng các kiểm định sự phù hợp và các tiêu chuẩn thông tin như AIC và BIC là rất quan trọng. Đầu tiên, chúng ta cần ước lượng các tham số của từng mô hình copula dựa trên dữ liệu phân tích. Sau đó, chúng ta sử dụng các kiểm định sự phù hợp như Anderson – Darling (A – D), Cramer – von Mises (Cv – M), hoặc Kolmogorov – Smornov (K – S) để đánh giá mức độ phù hợp của từng mô hình với dữ liệu thực tế. Các mô hình copula với giá trị p cao nhất hoặc không bị bác bỏ giả thuyết phù hợp được xem là phù hợp tốt hơn.

Đối với tiêu chuẩn thông tin, AIC và BIC đánh giá sự phù hợp của mô hình dựa trên hiệu suất phù hợp với dữ liệu và số lượng tham số. AIC đo lường mức độ mất mát thông tin trong mô hình và ưu tiên các mô hình có giá trị AIC thấp nhất. BIC, một biến thể của AIC, cố gắng tránh overfitting bằng cách ưa thích các mô hình đơn giản hơn và tìm kiếm mô hình có giá trị BIC thấp nhất, điều này được điều chỉnh dựa trên kích thước mẫu.

Tiếp cận này kết hợp cả kiểm định sự phù hợp thống kê và tiêu chuẩn thông tin để lựa chọn mô hình copula phù hợp nhất với dữ liệu. Điều này đảm bảo rằng mô hình copula được lựa chọn không bị overfitting hoặc underfitting và phù hợp tốt với dữ liệu thực tế.

3.3 Thiết kế nghiên cứu

Với mục đích nghiên cứu mức độ phản ứng của TTCK Việt Nam đối với sự biến động của TTCK Nhật Bản qua các cú sốc kinh tế, em quyết định thiết kế bài nghiên cứu như sau:

• Phần 1: Ước lượng mô hình Copula tìm ra cấu trúc phụ thuộc của TTCK Việt Nam và TTCK Nhật Bản trên toàn bộ dữ liệu với khoảng thời gian từ đầu 2018 đến đầu tháng 7/2024 nhằm đánh giá tổng quan mức độ phản ứng của TTCK Việt Nam xuyên suốt các thời kỳ.

• Phần 2: Nhằm so sánh sự khác biệt giữa mức độ phụ thuộc qua các cú sốc kinh tế và qua giai đoạn ổn định, em tiến hành ước lượng Copula cho từng giai đoạn thời gian cụ thể như đã chia ở phần Dữ liệu nghiên cứu..

Chương 4: Kết quả nghiên cứu

Trước khi đi vào nghiên cứu em tiến hành xử lý dữ liệu thô và tính toán tỷ suất lợi nhuận cho từng biến

library(openxlsx) 
library(quantmod)
library(tidyverse)
library(xts)
library(PerformanceAnalytics)
library(forecast)
library(rugarch)
library(goftest)
library(VineCopula)
library(tseries)
library(FinTS)
library(ggplot2)
library(corrplot)
#Nhập dữ liệu
a1 <- read.xlsx('/Users/phamxuanhoan/Downloads/Dữ liệu mô hình1 (1).xlsx')
a2 <- read.xlsx('/Users/phamxuanhoan/Downloads/Dữ liệu mô hình Nikkei (1).xlsx')
a1$Ngày <- as.Date(a1$Ngày)
a2$Ngày <- as.Date(a2$Ngày)
class(a1$Ngày)

## [1] "Date"

class(a2$Ngày)

## [1] "Date"

##Khai báo dữ liệu chuỗi thời gian, sắp xếp theo thời gian
VNI <- xts(a1[,-1], order.by = a1$Ngày)
NI225 <- xts(a2[,-1],order.by = a2$Ngày)
#Gộp hai chuỗi giá dựa trên những ngày trùng nhau
dt <- merge.xts(NI225,VNI,join = 'inner')

#Tính tỷ suất lợi nhuận theo công thức logt - logt-1
rt <- CalculateReturns(dt,method = 'log')
rt <- rt[-1]

4.1 Thống kê mô tả

Thống kê mô tả

library(kableExtra)
library(moments)
rt_df <- as.data.frame(rt)
a <- rt_df %>% summarise(Min = min(VNI),
                      Max = max(VNI),
                      Mean = mean(VNI),
                      StDev = sd(VNI),
                      Skewness = skewness(VNI),
                      Kurtosis = kurtosis(VNI))
b <- rt_df %>% summarise(Min = min(NI225),
                      Max = max(NI225),
                      Mean = mean(NI225),
                      StDev = sd(NI225),
                      Skewness = skewness(NI225),
                      Kurtosis = kurtosis(NI225))
m <- rbind(a,b)
rownames(m) <- c('VNIndex','Nikkei')
kable(m, format = 'pandoc', caption = 'Bảng 2: Thống kê mô tả chuỗi TSLN', table.attr = "style='width:100%;'") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bảng 2: Thống kê mô tả chuỗi TSLN

	Min	Max	Mean	StDev	Skewness	Kurtosis
VNIndex	-0.0852948	0.1217715	0.0001013	0.0139477	-0.4319986	11.576293
Nikkei	-0.0627357	0.0773137	0.0003820	0.0126139	-0.0440969	6.492025

Bảng 2 cho thấy kết quả thống kê mô tả của hai chuỗi tỷ suất sinh lợi của chỉ số VNINDEX đại diện cho TTCK Việt Nam và chỉ số NI225 đại diện cho TTCK Nhật Bản. Kết quả cho thấy, giá trị trung bình tỷ suất lợi nhuận của VNINDEX là 0.0001013 và của NI225 là -0.0003820. Chỉ số VNI có độ biến động (0.0139477) lớn hơn chỉ số NI225 (0.0126139), điều này chứng tỏ TTCK Việt Nam dao động khá lớn và không ổn định so với TTCK Nhật Bản, qua đó cho thấy rủi ro của TTCK Việt Nam cao hơn rất nhiều. Hệ số Skewness của cả hai chỉ số đều bé hơn 0 chứng tỏ rằng phân phối của cả hai chuỗi lợi suất đều lệch trái và phân phối của chuỗi VNI có đuôi dài hơn. Bên cạnh đó, các giá trị Kurtosis của cả hai đều lớn hơn rất nhiều so với mức 3, điều này chỉ ra rằng phân phối của cả hai chỉ số có đuôi dài và nhiều giá trị cực trị hơn so với phân phối chuẩn. Ta có thể trực quan hóa các mô tả cơ bản của hai chỉ số lợi suất thông qua biểu đồ hộp (Box-plot).

pivot_longer(rt_df, cols = everything(), names_to = "TSLN", values_to = "Value") %>% ggplot(aes(x = TSLN, y = Value)) +
  geom_boxplot(fill = "#FF66CC", color = "#FF3366") +
  theme_minimal() +
  labs(title = "Hình 7: Biểu đồ hộp của TSLN VNI và NI225",
       x = "Chỉ số",
       y = "Tỷ suất sinh lợi")

Biểu đồ hộp trên cho thấy rằng cả hai chỉ số NI225 và VNI đều có tỷ suất sinh lợi chủ yếu tập trung quanh giá trị trung vị, với một số điểm ngoại lệ ở cả hai phía, điều này chỉ ra rằng có một số giá trị tỷ suất sinh lợi bất thường trong dữ liệu, nguyên nhân có lẽ đến từ những cú sốc kinh tế từ đại dịch Covid-19 hay cuộc xung đột Nga-Ukraine.

Các kiểm định

#Kiểm định tính dừng
adf_vni <- adf.test(rt$VNI)
adf_ni225 <- adf.test(rt$NI225)
#Kiểm định phân phối chuẩn
jq_vni <- jarque.bera.test(rt$VNI)
jq_ni225 <- jarque.bera.test(rt$NI225)
# Ước lượng và trích xuất phần dư từ mô hình ARMA tối ưu
arima_VNI <- autoarfima(rt$VNI,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

arima_NI225 <- autoarfima(rt$NI225,ar.max = 2, ma.max = 2, criterion = 'AIC', method = "full")

re_VNI <- arima_VNI$fit@fit$residuals
re_NI225 <- arima_NI225$fit@fit$residuals
#Kiểm định tương quan chuỗi bậc 2 cho phần dư
lj_vni <- Box.test(re_VNI,type = 'Ljung-Box', lag = 2)
lj_ni225 <- Box.test(re_NI225,type = 'Ljung-Box', lag = 2)

#Kiểm định tương quan chuỗi bậc 2 cho phần dư bình phương
lj_vni2 <- Box.test(re_VNI^2,type = 'Ljung-Box', lag = 2)
lj_ni2252 <- Box.test(re_NI225^2,type = 'Ljung-Box', lag = 2)

#Kiểm định hiệu ứng ARCH
ar_vni <- ArchTest(re_VNI, lags = 2)
ar_ni225 <- ArchTest(re_NI225, lags = 2)

#Trình bày kết quả
test_result <- data.frame(Test = c("ADF","J-B","Q(2)","Q(2)^2", "ARCH(2)"),
                          P_value_VNI = c(adf_vni$p.value, jq_vni$p.value, lj_vni$p.value, lj_vni2$p.value, ar_vni$p.value),
                          P_value_NI225 = c(adf_ni225$p.value, jq_ni225$p.value, lj_ni225$p.value, lj_ni2252$p.value, ar_ni225$p.value))
kable(test_result, 
  caption = "Bảng 3: Kết quả các kiểm định", 
  label = 'Ghi chú: Q (2) and Q2 (2) lần lượt là kiểm định Ljung-Box Q2 cho tương quan chuỗi bậc 2 của phần dư và bình phương phần dư của lợi suất', 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 3: Kết quả các kiểm định

Test	P_value_VNI	P_value_NI225
ADF	0.0100000	0.0100000
J-B	0.0000000	0.0000000
Q(2)	0.6837059	0.9882016
Q(2)^2	0.0000000	0.0000000
ARCH(2)	0.0000000	0.0000000

Bảng 3 trình bày các kết quả kiểm định của chuỗi tỷ suất sinh lợi và phần dư mô hình ARMA của nó. Kết quả cho thấy cả hai chuỗi tỷ suất đều là chuỗi dừng và không có phân phối chuẩn dựa trên kiểm định Augmented – Dickey – Fuller (ADF) và Jarque – Bera. Kiểm định Ljung-Box chứng tỏ phần dư của chuỗi VNI không có hiện tượng tương quan chuỗi ở cả chuỗi gốc và bình phương. Với phần dư chuỗi NI225, kiểm định Ljung-Box cho chuỗi gốc (Q(2)) cho thấy không có hiện tượng tương quan chuỗi và kiểm định cho chuỗi bình phương (Q(2)^2) cho thấy hiện tượng tương quan chuỗi. Hiệu ứng ARCH không tồn tại ở phần dư của chuỗi NI225 và chuỗi VNI.

4.2 Hệ số tương quan

Trong phạm vi bài tiểu luận này, em thực hiện tính toán các hệ số tương quan bao gồm tương quan tuyến tính Pearson, tương quan hạng Spearman và Kendall.

pearson <- cor(rt$VNI,rt$NI225, method="pearson")
spearman <- cor(rt$VNI,rt$NI225, method="spearman")
kendall <- cor(rt$VNI,rt$NI225, method="kendall")

#Trình bày kết quả
relat <- data.frame('Tương quan' = 'VNI-NI225',
                    Pearson = pearson[1,1],
                    spearman = spearman[1,1],
                    Kendall = kendall[1,1])
kable(relat, 
      col.names = c("Phương pháp", "Pearson", "Spearman","Kendall"),
      caption = "Bảng 4: Kết quả hệ số tương quan",
      format = 'pandoc',
      align = c("l", "c", "c")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Bảng 4: Kết quả hệ số tương quan

Phương pháp	Pearson	Spearman	Kendall
VNI-NI225	0.294234	0.2772344	0.1910552

par(mfrow = c(1,2))

corr <- cor(rt)
library(ggcorrplot)
ggcorrplot(corr, hc.order = TRUE,
   outline.col = "white",
   ggtheme = ggplot2::theme_gray,
   colors = c("#6D9EC1", "white", "#E46726"),
   lab = TRUE, lab_col = 'white',title = 'Hình 8: Trực quan hóa hệ số tương quan với phương pháp Pearson')

chart.Correlation(rt, histogram=TRUE, pch=19)

Kết quả tính toán trên cho thấy hệ số tương quan của 3 phương pháp có mức chênh lệch nhau khá lớn, hệ số tương quan tuyến tính Pearson cho thấy mức độ tương quan lớn nhất với giá trị 0.29. Lưu ý rằng, hệ số tương quan tuyến tính Person giả định rằng các chuỗi lợi suất có phân phối chuẩn, nhưng thực tế thì lại không. Do đó, việc dùng hệ số tương quan để ước tính mối tương quan giữa 2 thị trường còn mang lại nhiều tranh cãi, và không thể giải thích được sự tương quan giữa các thị trường khi cả 2 biến động cực biên. Ngoài ra, hệ số tương quan hạng Spearman và Kendall không yêu cầu các chuỗi lợi suất có phân phối chuẩn. Hai hệ số tương quan này chưa phản ánh thông tin hoặc các cú sốc gây biến động thị trường, và chưa xem xét đến sự phụ thuộc đuôi của các thị trường. Do đó, để đánh giá khách quan hơn về mối quan hệ giữa hai thị trường, chúng ta sẽ tiến hành phương pháp Copula.

4.3 Phân tích biến động chuỗi tỷ suất sinh lợi

par(mfrow = c(1,2))
a11 <- CalculateReturns(a1, method = 'log')

ggplot(a1, aes(x = Ngày, y= a11$VNIndex))+
  geom_line()+
  labs(title = "Hình 9: Biến động theo ngày của TSLN VNINDEX",x = "Ngày", y="Tỷ suất sinh lợi")+
  theme(plot.title = element_text(hjust = 0.5))

a22 <- CalculateReturns(a2, method = 'log')

ggplot(a2, aes(x = Ngày, y= a22$Nikkei))+
  geom_line()+
  labs(title = "Hình 10: Biến động theo ngày của TSLN NI225",x = "Ngày", y="Tỷ suất sinh lợi")+
  theme(plot.title = element_text(hjust = 0.5))

Biến động của chuỗi tỷ suất sinh lợi ở chỉ số VNINDEX và chỉ số NI225 trong giai đoạn 2018 - 7/2024 được trình bày ở bảng trên. Nhìn chung, chỉ số NI225 cho thấy mức độ dao động mạnh hơn chỉ số VNI, điều này chứng tỏ rằng rủi ro ở TTCK Nhật Bản là lớn hơn rất nhiều so với TTCK Việt Nam. Lý do của sự biến động mạnh này là trong khoảng thời gian nghiên cứu, có xảy ra các cú sốc kinh tế do đại dịch Covid và cuộc xung đột Nga - Ukraine gây ra. Ở chỉ số NI225, tỷ suất sinh lợi dao động mạnh nhất nằm trong khoảng thời gian từ năm 2021 đến giữa năm 2022 khi lợi nhuận giảm tối đa ở mức 0.05, đây cũng là khoảng thời gian mà dịch Covid-19 và cuộc xung đột diễn ra. Các giai đoạn sau hai cú sốc cho thấy mức dao động dần ổn định hơn và ít các biến động cực biên. Với chỉ số VNI, mức dao động của tỷ suất sinh lợi trong khoảng thời gian nghiên cứu là khá ổn định, chỉ tồn tại các biến động cực biên ở khoảng thời gian đầu năm 2022, khi mà cuộc xung đột Nga - Ukraine vừa nổ ra. Điều này cho thấy có tồn tại sự ảnh hưởng của 2 cú sốc kinh tế đến thị trường chứng khoán ở cả Nhật Bản và Việt Nam.

Để có thể dễ dàng xác định rõ xu hướng biến động cùng với nhau giữa hai TTCK, chúng ta sẽ sử dụng biểu đồ Scatter với đường hồi quy để biểu diễn mối liên hệ tương quan giữa hai TTCK một cách trực quan

rt %>% ggplot(aes(VNI,NI225)) +
              geom_point(col = '#00CC66',shape = TRUE) + 
              geom_smooth(method = 'lm',se = T, col = '#CCFFFF') + #thêm đường hồi quy với phương pháp hồi quy tuyến tính  
              labs(title = 'Hình 11: Biểu đồ Scatter TSLN của NI225 và VNI', x = 'VNI', y = 'NI225')

Biểu đồ Scatter trên cho thấy hầu hết số lượng điểm dữ liệu nằm tập trung gần vị trí 0, đường hồi quy dốc lên chứng tỏ mối tương quan dương giữa tỷ suất sinh lợi của VNINDEX và NI225. Các điểm dữ liệu nằm rải rác xung quanh đường xu hướng, cho thấy rằng mặc dù có xu hướng tăng chung, nhưng vẫn có sự biến động lớn giữa các giá trị cụ thể của VNI và NI225.

4.4 Kết quả ước lượng cấu trúc phụ thuộc dựa trên phương pháp Copula có điều kiện

4.4.1 Xác định mô hình phân phối biên

Để ước lượng các tham số của mô hình copula, trước tiên cần xác định mô hình phân phối biên phù hợp nhất để mô tả chính xác đặc điểm biến động của từng chuỗi lợi suất chứng khoán. Các mô hình phân phối biên của sai số chuẩn hóa của mỗi chuỗi lợi suất là yếu tố đầu vào của copula hai biến. Mô hình biên tổng quát của các chuỗi lợi suất có dạng ARMA(p,q)-GJR-GARCH(r,m). Độ trễ p và q trong mô hình ARMA(p,q) được xác định dựa trên giá trị nhỏ nhất của tiêu chuẩn thông tin AIC, ước tính thông qua phương pháp ước lượng hợp lý cực đại (MLE). Đây là mô hình xác định độ trễ và bước nhảy tối ưu cho chuỗi dữ liệu. Hiệu ứng ARCH được phát hiện ở tất cả các chuỗi lợi suất thông qua kiểm định Lagrange Multiplier (LM). Mô hình biên ARMA(p,q)-GJR-GARCH(r,m) phù hợp nhất được chọn từ một tập hợp gồm 20 mô hình biên khác nhau, tương ứng với các kết hợp bậc của tham số r, m và độ trễ p, q tối đa 2; cùng với 5 dạng hàm phân phối biên của chuỗi phần dư cho mỗi chuỗi lợi suất. Các dạng hàm phân phối này bao gồm: phân phối chuẩn (Normal), Student-t, Skewed student-t, phân phối lỗi tổng quát (GED) và phân phối lỗi tổng quát lệch (sGED). Mô hình phân phối biên phù hợp nhất cho mỗi chuỗi lợi suất được chọn dựa trên giá trị nhỏ nhất của tiêu chuẩn thông tin AIC, BIC, SIC và HQIC.

#Setup garch(1,1) and estimate for VNI coef
#Norm
spec_VNI_ugarch11_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch11_norm <- ugarchfit(spec = spec_VNI_ugarch11_norm, data = rt$VNI)

#Std
spec_VNI_ugarch11_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch11_std <- ugarchfit(spec = spec_VNI_ugarch11_std, data = rt$VNI)

#Sstd
spec_VNI_ugarch11_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch11_sstd <- ugarchfit(spec = spec_VNI_ugarch11_sstd, data = rt$VNI)

#ged
spec_VNI_ugarch11_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch11_ged <- ugarchfit(spec = spec_VNI_ugarch11_ged, data = rt$VNI)

#Sged
spec_VNI_ugarch11_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch11_sged <- ugarchfit(spec = spec_VNI_ugarch11_sged, data = rt$VNI)

#Setup garch(1,2) and estimate for VNI coef
#norm
spec_VNI_ugarch12_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch12_norm <- ugarchfit(spec = spec_VNI_ugarch12_norm, data = rt$VNI)

#std
spec_VNI_ugarch12_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch12_std <- ugarchfit(spec = spec_VNI_ugarch12_std, data = rt$VNI)

#sstd
spec_VNI_ugarch12_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch12_sstd <- ugarchfit(spec = spec_VNI_ugarch12_sstd, data = rt$VNI)

#ged
spec_VNI_ugarch12_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch12_ged <- ugarchfit(spec = spec_VNI_ugarch12_ged, data = rt$VNI)

#sged
spec_VNI_ugarch12_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch12_sged <- ugarchfit(spec = spec_VNI_ugarch12_sged, data = rt$VNI)

#Setup garch(2,1) and estimate for VNI coef
#norm
spec_VNI_ugarch21_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch21_norm <- ugarchfit(spec = spec_VNI_ugarch21_norm, data = rt$VNI)

#std
spec_VNI_ugarch21_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch21_std <- ugarchfit(spec = spec_VNI_ugarch21_std, data = rt$VNI)

#sstd
spec_VNI_ugarch21_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch21_sstd <- ugarchfit(spec = spec_VNI_ugarch21_sstd, data = rt$VNI)

#ged
spec_VNI_ugarch21_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch21_ged <- ugarchfit(spec = spec_VNI_ugarch21_ged, data = rt$VNI)

#sged
spec_VNI_ugarch21_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch21_sged <- ugarchfit(spec = spec_VNI_ugarch21_sged, data = rt$VNI)

#Setup garch(2,2) and estimate for VNI coef
#norm
spec_VNI_ugarch22_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
  distribution.model = "norm")

VNI_garch22_norm <- ugarchfit(spec = spec_VNI_ugarch22_norm, data = rt$VNI)

#std
spec_VNI_ugarch22_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "std")

VNI_garch22_std <- ugarchfit(spec = spec_VNI_ugarch22_std, data = rt$VNI)

#sstd
spec_VNI_ugarch22_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sstd")

VNI_garch22_sstd <- ugarchfit(spec = spec_VNI_ugarch22_sstd, data = rt$VNI)

#ged
spec_VNI_ugarch22_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "ged")

VNI_garch22_ged <- ugarchfit(spec = spec_VNI_ugarch22_ged, data = rt$VNI)

#sged
spec_VNI_ugarch22_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
  distribution.model = "sged")

VNI_garch22_sged <- ugarchfit(spec = spec_VNI_ugarch22_sged, data = rt$VNI)

#Setup garch(1,1) and estimate for NI225 coef
#norm
spec_NI225_ugarch11_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

NI225_garch11_norm <- ugarchfit(spec = spec_NI225_ugarch11_norm, data = rt$NI225)

#std
spec_NI225_ugarch11_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

NI225_garch11_std <- ugarchfit(spec = spec_NI225_ugarch11_std, data = rt$NI225)

#sstd
spec_NI225_ugarch11_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

NI225_garch11_sstd <- ugarchfit(spec = spec_NI225_ugarch11_sstd, data = rt$NI225)

#ged
spec_NI225_ugarch11_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

NI225_garch11_ged <- ugarchfit(spec = spec_NI225_ugarch11_ged, data = rt$NI225)

#sged
spec_NI225_ugarch11_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

NI225_garch11_sged <- ugarchfit(spec = spec_NI225_ugarch11_sged, data = rt$NI225)

#Setup garch(1,2) and estimate for NI225 coef
#norm
spec_NI225_ugarch12_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

NI225_garch12_norm <- ugarchfit(spec = spec_NI225_ugarch12_norm, data = rt$NI225)

#std
spec_NI225_ugarch12_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

NI225_garch12_std <- ugarchfit(spec = spec_NI225_ugarch12_std, data = rt$NI225)

#sstd
spec_NI225_ugarch12_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

NI225_garch12_sstd <- ugarchfit(spec = spec_NI225_ugarch12_sstd, data = rt$NI225)

#ged
spec_NI225_ugarch12_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

NI225_garch12_ged <- ugarchfit(spec = spec_NI225_ugarch12_ged, data = rt$NI225)

#sged
spec_NI225_ugarch12_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(1, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

NI225_garch12_sged <- ugarchfit(spec = spec_NI225_ugarch12_sged, data = rt$NI225)

#Setup garch(2,1) and estimate for NI225 coef
#norm
spec_NI225_ugarch21_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

NI225_garch21_norm <- ugarchfit(spec = spec_NI225_ugarch21_norm, data = rt$NI225)

#std
spec_NI225_ugarch21_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

NI225_garch21_std <- ugarchfit(spec = spec_NI225_ugarch21_std, data = rt$NI225)

#sstd
spec_NI225_ugarch21_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

NI225_garch21_sstd <- ugarchfit(spec = spec_NI225_ugarch21_sstd, data = rt$NI225)

#ged
spec_NI225_ugarch21_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

NI225_garch21_ged <- ugarchfit(spec = spec_NI225_ugarch21_ged, data = rt$NI225)

#sged
spec_NI225_ugarch21_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

NI225_garch21_sged <- ugarchfit(spec = spec_NI225_ugarch21_sged, data = rt$NI225)

#Setup garch(2,2) and estimate for NI225 coef
#norm
spec_NI225_ugarch22_norm <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "norm")

NI225_garch22_norm <- ugarchfit(spec = spec_NI225_ugarch22_norm, data = rt$NI225)

#std
spec_NI225_ugarch22_std <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "std")

NI225_garch22_std <- ugarchfit(spec = spec_NI225_ugarch22_std, data = rt$NI225)

#sstd
spec_NI225_ugarch22_sstd <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sstd")

NI225_garch22_sstd <- ugarchfit(spec = spec_NI225_ugarch22_sstd, data = rt$NI225)

#ged
spec_NI225_ugarch22_ged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "ged")

NI225_garch22_ged <- ugarchfit(spec = spec_NI225_ugarch22_ged, data = rt$NI225)

#sged
spec_NI225_ugarch22_sged <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = c(2, 2)),
  mean.model = list(armaOrder = c(0, 2), include.mean = TRUE),
  distribution.model = "sged")

NI225_garch22_sged <- ugarchfit(spec = spec_NI225_ugarch22_sged, data = rt$NI225)

#Choose optimization model for vni coef
VNI.model.list <- list(garch11n = VNI_garch11_norm, garch11t = VNI_garch11_std, garch11st = VNI_garch11_sstd, 
                       garch11g = VNI_garch11_ged, garch11sg = VNI_garch11_sged, 
                       garch12n = VNI_garch12_norm, garch12t = VNI_garch12_std, garch12st = VNI_garch12_sstd, 
                       garch12g = VNI_garch12_ged, garch12sg = VNI_garch12_sged,
                       garch21n = VNI_garch21_norm, garch21t = VNI_garch21_std, garch21st = VNI_garch21_sstd, 
                       garch21g = VNI_garch21_ged, garch21sg = VNI_garch21_sged,
                       garch22n = VNI_garch22_norm, garch22t = VNI_garch22_std, garch22st = VNI_garch22_sstd, 
                       garch22g = VNI_garch22_ged, garch22sg = VNI_garch22_sged
                       ) 

VNI.info.mat <- sapply(VNI.model.list, infocriteria)

rownames(VNI.info.mat) <- rownames(infocriteria(VNI_garch11_norm))

VNI.inds <- which(VNI.info.mat == min(VNI.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.VNI <- colnames(VNI.info.mat)[VNI.inds[,2]]


#Choose optimization model for NI225 coef
NI225.model.list <- list(garch11n = NI225_garch11_norm, garch11t = NI225_garch11_std, garch11st = NI225_garch11_sstd, 
                       garch11g = NI225_garch11_ged, garch11sg = NI225_garch11_sged, 
                       garch12n = NI225_garch12_norm, garch12t = NI225_garch12_std, garch12st = NI225_garch12_sstd, 
                       garch12g = NI225_garch12_ged, garch12sg = NI225_garch12_sged,
                       garch21n = NI225_garch21_norm, garch21t = NI225_garch21_std, garch21st = NI225_garch21_sstd, 
                       garch21g = NI225_garch21_ged, garch21sg = NI225_garch21_sged,
                       garch22n = NI225_garch22_norm, garch22t = NI225_garch22_std, garch22st = NI225_garch22_sstd, 
                       garch22g = NI225_garch22_ged, garch22sg = NI225_garch22_sged
                       )

NI225.info.mat <- sapply(NI225.model.list, infocriteria)

rownames(NI225.info.mat) <- rownames(infocriteria(NI225_garch11_norm))

NI225.inds <- which(NI225.info.mat == min(NI225.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.NI225 <- colnames(NI225.info.mat)[NI225.inds[,2]]

#Trình bày kết quả
mar_model <- data.frame(
  rate = c("VNI", "NI225"),
  Format = c("ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t", "ARMA(0,2)-GJR-GARCH(1,1)-Student t"))

# Render the table
kable(mar_model, col.names = c("Tỷ suất sinh lợi", "Dạng mô hình phân phối biên"), 
      caption = "Bảng 5: Mô hình phân phối biên tối ưu", format = "pandoc", 
      table.attr = "style='width:100%;'") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bảng 5: Mô hình phân phối biên tối ưu

Tỷ suất sinh lợi	Dạng mô hình phân phối biên
VNI	ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t
NI225	ARMA(0,2)-GJR-GARCH(1,1)-Student t

Dựa trên tiêu chuẩn giá trị AIC nhỏ nhất, kết quả thu được dạng mô hình phân phối biên cho tỷ suất sinh lợi của VNI và NI225 lần lượt là mô hình ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student t và ARMA(0,2)-GJR-GARCH(1,1)-Student t. Điều này một lần nữa xác nhận rằng, hàm phân phối biên của tất cả các chuỗi lợi suất chứng khoán được xem xét không tuân theo phân phối chuẩn.

#Trình bày kết quả

extract_garch_results <- function(fit) {
  coef <- coef(fit)
  se <- fit@fit$se.coef
  pvalues <- 2 * (1 - pnorm(abs(fit@fit$tval))) # tính giá trị p từ giá trị t
  results <- cbind(coef, se, pvalues)
  colnames(results) <- c("Estimate", "Std. Error", "Pr(>|z|)")
  return(results)
}

fit1 <- extract_garch_results(VNI_garch11_sstd)
fit2 <- extract_garch_results(NI225_garch11_std)

kable(as.data.frame(fit1), 
  caption = "Bảng 5: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Skewed Student của biến VNI", 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 5: Kết quả mô hình ARMA(1,1)-GJR-Garch(1,1)-Skewed Student của biến VNI

	Estimate	Std. Error	Pr(>|z|)
mu	0.0004476	0.0002877	0.1197833
ar1	-0.9347537	0.0269390	0.0000000
ma1	0.9470713	0.0240976	0.0000000
omega	0.0000194	0.0000058	0.0008854
alpha1	0.0136137	0.0235398	0.5630424
beta1	0.7324710	0.0537623	0.0000000
gamma1	0.3379362	0.0950513	0.0003775
skew	0.8700608	0.0321430	0.0000000
shape	3.4294383	0.3419997	0.0000000

Bảng 5 kết quả mô hình ARMA(1,1)-GJR-GARCH(1,1)-Skewed Student cho tỷ suất sinh lợi của chỉ số VNI cho thấy một số đặc điểm quan trọng về cấu trúc và động thái của dữ liệu. Tham số AR(1) có giá trị âm gần 1 và có ý nghĩa thống kê thấp, cho thấy sự tự duy trì mạnh mẽ của lợi nhuận qua các thời kỳ. Ngược lại, tham số MA(1) có giá trị âm và cũng rất có ý nghĩa, điều này ngụ ý rằng các cú sốc quá khứ không ảnh hưởng tiêu cực mạnh mẽ đến lợi nhuận hiện tại. Trong phần mô hình GARCH, tham số omega rất nhỏ nhưng có ý nghĩa, biểu thị mức độ biến động cơ bản thấp. Tham số beta1 cao và có ý nghĩa, cho thấy sự duy trì mạnh mẽ của biến động qua thời gian. Tham số alpha1 không có ý nghĩa thống kê, nghĩa là các cú sốc quá khứ không ảnh hưởng nhiều đến biến động hiện tại. Tuy nhiên, tham số gamma1 lại có ý nghĩa, cho thấy hiệu ứng đòn bẩy, nơi mà các cú sốc tiêu cực có ảnh hưởng mạnh hơn đến biến động so với các cú sốc tích cực có cùng độ lớn. Cuối cùng, tham số skew và shape đều có ý nghĩa thống kê, chỉ ra rằng phân phối sai số có độ lệch và đuôi dày hơn so với phân phối chuẩn. Cụ thể, độ lệch dương cho thấy phân phối bị lệch về phía phải, trong khi giá trị shape cho thấy sự hiện diện của các đuôi nặng, ám chỉ rằng các giá trị cực đoan xảy ra thường xuyên hơn so với phân phối chuẩn.

kable(as.data.frame(fit2), 
  caption = "Bảng 6: Kết quả mô hình ARMA(0,2)-GJR-Garch(1,1)-Student của biến NI225", 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 6: Kết quả mô hình ARMA(0,2)-GJR-Garch(1,1)-Student của biến NI225

	Estimate	Std. Error	Pr(>|z|)
mu	0.0004555	0.0002970	0.1251750
ma1	-0.0116990	0.0271076	0.6660486
ma2	0.0673197	0.0261084	0.0099238
omega	0.0000134	0.0000003	0.0000000
alpha1	0.0000001	0.0081318	0.9999881
beta1	0.8236105	0.0163566	0.0000000
gamma1	0.1872201	0.0376615	0.0000007
shape	5.7915527	0.9117449	0.0000000

Với mô hình ARMA(0,2)-GJR-GARCH(1,1)-Student cho tỷ suất sinh lợi của chỉ số NI225 cho thấy tham số MA(1) không có ý nghĩa thống kê và MA(2) có ý nghĩa thống kê, nghĩa là các cú sốc quá khứ gần đây một phần có ảnh hưởng đáng kể đến lợi nhuận hiện tại. Đối với mô hình GARCH, tham số omega rất nhỏ nhưng có ý nghĩa thống kê, biểu thị mức độ biến động cơ bản thấp. Tương tự như trên, tham số alpha1 của mô hình NI225 cũng không có ý nghĩa thống kê, cho thấy các cú sốc quá khứ không ảnh hưởng nhiều đến biến động hiện tại. Ngược lại, tham số beta1 cao và có ý nghĩa, cho thấy sự duy trì mạnh mẽ của biến động qua thời gian. Tham số gamma1 có ý nghĩa thống kê, chỉ ra hiệu ứng đòn bẩy, nơi mà các cú sốc tiêu cực có ảnh hưởng mạnh hơn đến biến động so với các cú sốc tích cực có cùng độ lớn. Bên cạnh đó, tham số shape cũng chỉ ra rằng phân phối sai số có đuôi dày hơn so với phân phối chuẩn, ám chỉ rằng các giá trị cực đoan xảy ra thường xuyên hơn. Tóm lại, kết quả mô hình cho thấy rằng mặc dù các thành phần trung bình và cú sốc quá khứ gần đây không có ảnh hưởng đáng kể đến lợi nhuận hiện tại, nhưng mức độ biến động có tính duy trì cao và có hiệu ứng đòn bẩy đáng kể. Phân phối sai số cũng cho thấy đuôi dày, ngụ ý rằng các biến động lớn xảy ra thường xuyên hơn so với phân phối chuẩn.

4.4.2 Kiểm định tính phù hợp của mô hình phân phối biên

Sau khi chọn mô hình phân phối biên phù hợp nhất, chúng ta tiến hành trích xuất phần dư chuẩn hóa \(z_{1t}|\Omega_{t-i}, z_{2t}|\Omega_{t-i}\) cho từng chuỗi lợi suất. Tiếp theo, chuỗi dữ liệu này được chuyển đổi thành các giá trị tích phân xác suất thông qua hàm phân phối thực nghiệm F, được đề xuất bởi Genest và cộng sự (1995). Để thỏa mãn điều kiện đầu vào của hàm copula, các biến \(u_t\) và \(v_t\) được giả định là độc lập và có phân phối đồng nhất trên [0,1]. Bảng 7 trình bày kết quả của các kiểm định Anderson-Darling, Cramer-von Mises (Cv-M), và Kolmogorov-Smirnov (K-S).

# make probability integration for residual of vni
VNI.res <- residuals(VNI_garch11_sstd)/sigma(VNI_garch11_sstd) #standarded resid

a <- fitdist(distribution = "sstd", VNI.res, control = list()) #fit distribution for this resid

u <- pdist(distribution = "sstd", q = VNI.res, mu = a$pars['mu'] , sigma = a$pars['sigma'],
           skew = a$pars['skew'], shape = a$pars['shape']) # make probability integration

# make probability integration for residual of NI225
NI225.res <- residuals(NI225_garch11_std)/sigma(NI225_garch11_std) #standarđ resid

b <- fitdist(distribution = "std", NI225.res, control = list()) #fit distribution for this resid

v <- pdist(distribution = "std", q = NI225.res, mu = b$pars['mu'] , sigma = b$pars['sigma'],
           skew = b$pars['skew'], shape = b$pars['shape']) # make probability integration


#test with Anderson-Darling
library(goftest)
ad_vni <- ad.test(u, "punif")
ad_NI225 <- ad.test(v, "punif")


#Kiểm định Cramer-von Mises
cvm_vni <- cvm.test(u, "punif")
cvm_NI225 <- cvm.test(v, "punif")


# Kiem dinh ks-test
ks_vni <- ks.test(u, "punif")
ks_NI225 <- ks.test(v, "punif")

#Trình bày kết quả
test <- data.frame(test = c('P_value VNI' ,'P_value NI225'),
                    AD = c(ad_vni$p.value, ad_NI225$p.value),
                    CVM = c(cvm_vni$p.value, cvm_NI225$p.value),
                   KS = c(ks_vni$p.value,ks_NI225$p.value))

kable(test, 
  caption = "Bảng 7: Kết quả các kiểm định cho mô hình phân phối biên", col.names = c("Các kiểm định","Anderson-Darling", "Cramer-von Mises","Kolmogorov-Smirnov"), 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 7: Kết quả các kiểm định cho mô hình phân phối biên

Các kiểm định	Anderson-Darling	Cramer-von Mises	Kolmogorov-Smirnov
P_value VNI	0.7823030	0.7073582	0.4427939
P_value NI225	0.6428772	0.6961288	0.8522831

Bảng 7 trình bày kết quả các kiểm định phân phối cho mô hình của hai biến VNI và NI225, với ba loại kiểm định: Anderson-Darling, Cramer-von Mises và Kolmogorov-Smirnov. Các giá trị P_value cao và lớn hơn mức ý nghĩa từ cả ba kiểm định cho thấy rằng không có bằng chứng thống kê để bác bỏ giả thuyết rằng dữ liệu của các chuỗi tỷ suất VNI và NI225 tuân theo phân phối được giả định trong mô hình. Cụ thể, giá trị P_value của kiểm định Anderson-Darling là 0.782303 cho VNI và 0.6428772 cho NI225, giá trị P_value của kiểm định Cramer-von Mises là 0.6073582 cho VNI và 0.6961288 cho NI225, và giá trị P_value của kiểm định Kolmogorov-Smirnov là 0.4427939 cho VNI và 0.8522831 cho NI225. Những kết quả này củng cố kết luận rằng mô hình Copula sẽ là phù hợp với dữ liệu của hai biến này.

4.4.3 Ước lượng tham số mô hình copula

Sau khi lựa chọn được mô hình biên phù hợp, em sẽ ước lượng mô hình Copula dựa trên phương pháp MLE, mô hình Copula tối ưu được lựa chọn dựa trên chỉ số AIC và BIC nhỏ nhất với hàm BiCopSelect nằm trong package VineCopula.

Sự phụ thuộc giữa TTCK Việt Nam đối với TTCK Nhật Bản trong điều kiện biến động bình thường

Copula họ Elip (Gauss và Studentt) cho phép mô tả mối quan hệ giữa hai thị trường trong điều kiện thị trường biến động bình thường. Bảng 8 trình bày kết quả tham số ước lượng của các hàm copula họ Elip là Gauss và Student và hệ số phụ thuộc đuôi thể hiện cấu trúc phụ thuộc giữa TTCK Việt Nam và TTCK Nhật Bản.

#Ước lượng
gau <- BiCopEst(u,v,family = 1, method = "mle", se = T, max.df = 10)
stu <- BiCopEst(u,v,family = 2, method = "mle", se = T, max.df = 10)
#Trình bày kết quả

est <- data.frame(mqh = c('VNI-NI225' ,'VNI-NI225'),
                    Copula = c('Gauss', 'Student'),
                    Thamso = c(gau$par, stu$par),
                  Thamso2 = c(gau$par2, stu$par2),
                   duoi = c(gau$taildep$lower,stu$taildep$lower))

kable(est, 
  caption = "Bảng 8: Kết quả ước lượng Copula họ Elip", col.names = c("Chỉ số","Copula", "Par","Par2", "Phụ thuộc đuôi"), 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 8: Kết quả ước lượng Copula họ Elip

Chỉ số	Copula	Par	Par2	Phụ thuộc đuôi
VNI-NI225	Gauss	0.2872437	0.000000	0.000000
VNI-NI225	Student	0.2865164	9.285405	0.037423

Kết quả ước lượng copula họ Elip trong điều kiện thị trường bình thường cho mối quan hệ giữa hai biến VNI và NI225, bao gồm hai loại copula: Gaussian và Student. Kết quả từ copula Gaussian cho thấy một mức độ tương quan dương giữa hai biến với giá trị Par là 0.2872437, nhưng không thể mô hình hóa phụ thuộc đuôi, như thể hiện qua giá trị phụ thuộc đuôi là 0. Ngược lại, copula Student không chỉ cho thấy mức độ tương quan dương tương tự với giá trị Par là 0.2865164, mà còn có khả năng mô hình hóa các mối quan hệ phi tuyến và phụ thuộc đuôi. Tham số Par2 của copula Student là 9.249898, cho thấy độ tự do của phân phối Student t, và giá trị phụ thuộc đuôi là 0.037423, biểu thị mức độ phụ thuộc đuôi dương giữa VNI và NI225. Điều này ngụ ý rằng có một xác suất nhất định rằng hai biến sẽ cùng trải qua các biến động cực đoan theo cùng một hướng. Cụ thể, nếu TTCK Nhật Bản biến động tăng (giảm) giá thì khả năng TTCK Việt Nam tăng (giảm) theo là khoảng 3.74%. Mối tương quan dương này cũng được Le Van Thu(2020) chứng minh trong bài nghiên cứu của mình.

Dựa trên kết quả ước lượng các tham số Copula, em tiến hành mô phỏng dữ liệu dựa trên các tham số của Copula tương ứng và trực quan hóa mối quan hệ của TTCk Việt Nam và TTCK Nhật Bản với phối cảnh hàm phân phối xác suất (PDF) của nó.

set.seed(123)
#Mô phỏng copula gauss
cop_nor <- normalCopula(param = gau$par, dim = 2)

#Mô phỏng copula student
cop_std <- tCopula(param = stu$par, dim = 2, df = 1)

nor_per <- persp_plot(cop_nor, 'norp.png', '#FF6666')
std_per <- persp_plot(cop_std, 'stdp.png', '#CCFF33')

legend <- legendGrob(
  labels = c("Gauss", "Student"), pch = 15,
  gp = gpar(col = c('#FF6666', '#CCFF33'), fill = c('#FF6666', '#CCFF33'))
)

grid.arrange( nor_per, std_per, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 3)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 12: Phối cảnh PDF của Copula họ Elip được ước lượng", 
                 gp = gpar(fontsize = 15, font = 2))
             )

Sự phụ thuộc giữa TTCK Việt Nam đối với TTCK Nhật Bản trong điều kiện biến động cực biên

Khi thị trường biến động cực biên, cấu trúc phụ thuộc giữa TTCK Việt Nam và TTCK Nhật Bản được thể hiện qua mô hình Copula phù hợp nhất dựa trên chỉ số AIC.

#Lựa chọn và ước lượng
opt_cop <- BiCopSelect(u,v,selectioncrit = "AIC",method = "mle")
est_opt_cop <- BiCopEst(u,v,family = opt_cop$family,method = "mle",max.df = 30)

#Trình bày kết quả
est_opt_cop_dt <- data.frame(mqh = c('VNI-NI225'),
                    Copula = est_opt_cop$familyname,
                    Thamso = est_opt_cop$par,
                  Thamso2 = est_opt_cop$par2,
                   duoiduoi = est_opt_cop$taildep$lower,
                  duoitren = est_opt_cop$taildep$upper,
                  tau = est_opt_cop$tau,
                  aic = est_opt_cop$AIC,
                  bic = est_opt_cop$BIC)

kable(est_opt_cop_dt, 
  caption = "Bảng 9: Kết quả ước lượng mô hình Copula tối ưu", col.names = c("Chỉ số","Copula", "Par","Par2", "Đuôi dưới", "Đuôi trên","Hệ số Kendall", "AIC","BIC"), 
  format = 'pandoc') %>% kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 9: Kết quả ước lượng mô hình Copula tối ưu

Chỉ số	Copula	Par	Par2	Đuôi dưới	Đuôi trên	Hệ số Kendall	AIC	BIC
VNI-NI225	Survival BB1	0.1251005	1.152925	0.1756798	0.0081822	0.1837008	-130.4951	-120.0441

Bảng 9 trình bày kết quả ước lượng mô hình copula BB7 cho mối quan hệ giữa hai chỉ số tỷ suất lợi nhuận VNI và NI225 trong giai đoạn từ đầu năm 2018 đến tháng 7/2024. Tham số Par là 0.1251005 và Par2 là 1.152925 biểu thị mức độ phụ thuộc và hình dạng của sự phụ thuộc giữa hai chỉ số này. Mức độ phụ thuộc đuôi dưới là 0.1756798, cho thấy khả năng mà cả hai chỉ số cùng giảm cực đoan là 1.756%, điều này ngụ ý rằng trong giai đoạn nghiên cứu, có một mức độ nhất định mà các biến cố giảm lớn cùng xảy ra ở cả hai thị trường nhưng không quá cao. Mức độ phụ thuộc đuôi trên là 0.0081822, biểu thị khả năng cả hai chỉ số cùng tăng cực đoan là 9.97%, cho thấy có sự liên kết mạnh hơn khi thị trường tăng mạnh, điều này có thể được hiểu là trong các giai đoạn thị trường tăng mạnh, cả hai chỉ số VNI và NI225 có xu hướng cùng tăng. Hệ số Kendall là 0.1837008, chỉ ra một mối quan hệ tích cực nhưng không quá mạnh giữa hai chỉ số, cho thấy rằng khi một thị trường tăng hoặc giảm, thị trường kia cũng có xu hướng di chuyển theo cùng một hướng nhưng với mức độ liên kết không quá chặt chẽ. Điều này có ý nghĩa quan trọng trong việc quản lý rủi ro và đa dạng hóa danh mục đầu tư. Các chỉ số AIC và BIC lần lượt là -130.4951 và -120.0441, cho thấy mô hình copula BB7 có độ phù hợp tốt với dữ liệu trong giai đoạn này. Giá trị âm của AIC và BIC chỉ ra rằng mô hình đã mô tả tốt mối quan hệ phụ thuộc giữa VNI và NI225, làm rõ rằng copula BB7 là một lựa chọn tối ưu để mô hình hóa mối quan hệ này. Điều này có thể giúp các nhà đầu tư và nhà quản lý rủi ro hiểu rõ hơn về mối quan hệ giữa hai chỉ số thị trường lớn này và đưa ra các quyết định đầu tư và phòng ngừa rủi ro một cách hiệu quả hơn. Kết quả này là tương đồng với nghiên cứu của Le Van Thu (2020) về mối quan hệ giữa TTCK Việt Nam đối với sự biến động của TTCK Nhật Bản

So sánh mức độ phản ứng của TTCK Việt Nam đối với sự biến động TTCK Nhật Bản

Tiếp theo em sẽ vận dụng mô hình copula có điều kiện (ARMA – GJR – GARCH copula) để so sánh cấu trúc và mức độ phụ thuộc đuôi giữa TTCK Việt Nam và Nhật Bản giai đoạn trước cú sốc, đại dịch COVID – 19, chiến tranh Nga – Ukraine, giai đoạn sau cú sốc nhằm xác định liệu có tồn tại hay không sự lây lan biến động các cú sốc kinh tế của TTCK Nhật Bản đến TTCK Việt Nam. Căn cứ xác định mốc thời gian cho mỗi giai đoạn đã được lập luận ở trong phần dữ liệu nghiên cứu. Dạng mô hình phân phối biên mô tả sự biến động giá của từng chuỗi lợi suất trong 4 được trình bày trong bảng 10.

#Trích xuất dữ liệu
truoc <- rt["2088-01-01/2089-12-31"]
covid <- rt["2090-01-01/2091-12-31"]
xung_dot <- rt["2092-02-24/2093-02-21"]
sau <- rt["2093-03-01/2094-07-01"]

#Estimate arma
arma1vni <- autoarfima(truoc$VNI, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma1NI225 <- autoarfima(truoc$NI225, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma2vni <- autoarfima(covid$VNI, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma2NI225 <- autoarfima(covid$NI225, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma3vni <- autoarfima(xung_dot$VNI, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma3NI225 <- autoarfima(xung_dot$NI225, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma4vni <- autoarfima(sau$VNI, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")
arma4NI225 <- autoarfima(sau$NI225, ar.max = 2, ma.max = 2,criterion = 'AIC', method = "full")

arma1vni <- c(2,0)
arma1NI225 <- c(0,0)
arma2NI225 <- c(0,2)
arma2vni <- c(2,2)
arma3NI225 <- c(2,2)
arma3vni <- c(2,2)
arma4NI225 <- c(1,0)
arma4vni <- c(2,2)

#Setup garch and estimate for VNI coef

fit_gjrgarch <- function(data, amorder, garchoder, dis) {
  spec <- ugarchspec(
  variance.model = list(model = "gjrGARCH", garchOrder = garchoder),
  mean.model = list(armaOrder = amorder, include.mean = TRUE),
  distribution.model = dis)
  ugarchfit(spec = spec, data = data)
}
#Giai đoạn trước cho VNI
garch11_nor_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,1), "norm")
garch11_std_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,1), "std")
garch11_sstd_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,1), "sstd")
garch11_ged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,1), "ged")
garch11_sged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,1), "sged")
garch22_nor_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,2), "norm")
garch22_std_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,2), "std")
garch22_sstd_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,2), "sstd")
garch22_ged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,2), "ged")
garch22_sged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,2), "sged")
garch12_nor_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,2), "norm")
garch12_std_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,2), "std")
garch12_sstd_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,2), "sstd")
garch12_ged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,2), "ged")
garch12_sged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(1,2), "sged")
garch21_nor_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,1), "norm")
garch21_std_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,1), "std")
garch21_sstd_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,1), "sstd")
garch21_ged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,1), "ged")
garch21_sged_vni1 <- fit_gjrgarch(truoc$VNI, arma1vni, c(2,1), "sged")
#Giai đoạn Covid cho VNI
garch11_nor_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,1), "norm")
garch11_std_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,1), "std")
garch11_sstd_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,1), "sstd")

garch11_ged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,1), "ged")

garch11_sged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,1), "sged")
garch22_nor_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,2), "norm")
garch22_std_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,2), "std")
garch22_sstd_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,2), "sstd")
garch22_ged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,2), "ged")
garch22_sged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,2), "sged")
garch12_nor_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,2), "norm")
garch12_std_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,2), "std")
garch12_sstd_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,2), "sstd")
garch12_ged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,2), "ged")
garch12_sged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(1,2), "sged")
garch21_nor_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,1), "norm")
garch21_std_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,1), "std")
garch21_sstd_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,1), "sstd")
garch21_ged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,1), "ged")
garch21_sged_vni2 <- fit_gjrgarch(covid$VNI, arma2vni, c(2,1), "sged")
#Giai đoạn xung đột
garch11_nor_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,1), "norm")
garch11_std_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,1), "std")
garch11_sstd_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,1), "sstd")
garch11_ged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,1), "ged")
garch11_sged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,1), "sged")
garch22_nor_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,2), "norm")
garch22_std_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,2), "std")
garch22_sstd_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,2), "sstd")
garch22_ged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,2), "ged")
garch22_sged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,2), "sged")
garch12_nor_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,2), "norm")
garch12_std_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,2), "std")
garch12_sstd_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,2), "sstd")
garch12_ged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,2), "ged")
garch12_sged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(1,2), "sged")
garch21_nor_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,1), "norm")
garch21_std_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,1), "std")
garch21_sstd_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,1), "sstd")
garch21_ged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,1), "ged")
garch21_sged_vni3 <- fit_gjrgarch(xung_dot$VNI, arma3vni, c(2,1), "sged")
#Giai đoạn sau cú sốc 
garch11_nor_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,1), "norm")
garch11_std_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,1), "std")
garch11_sstd_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,1), "sstd")
garch11_ged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,1), "ged")
garch11_sged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,1), "sged")
garch22_nor_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,2), "norm")
garch22_std_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,2), "std")
garch22_sstd_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,2), "sstd")
garch22_ged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,2), "ged")
garch22_sged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,2), "sged")
garch12_nor_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,2), "norm")
garch12_std_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,2), "std")
garch12_sstd_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,2), "sstd")
garch12_ged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,2), "ged")
garch12_sged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(1,2), "sged")
garch21_nor_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,1), "norm")
garch21_std_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,1), "std")
garch21_sstd_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,1), "sstd")
garch21_ged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,1), "ged")
garch21_sged_vni4 <- fit_gjrgarch(sau$VNI, arma4vni, c(2,1), "sged")
#Cho NI225
garch11_nor_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,1), "norm")
garch11_std_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,1), "std")
garch11_sstd_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,1), "sstd")
garch11_ged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,1), "ged")
garch11_sged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,1), "sged")
garch22_nor_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,2), "norm")
garch22_std_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,2), "std")
garch22_sstd_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,2), "sstd")
garch22_ged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,2), "ged")
garch22_sged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,2), "sged")
garch12_nor_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,2), "norm")
garch12_std_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,2), "std")
garch12_sstd_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,2), "sstd")
garch12_ged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,2), "ged")
garch12_sged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(1,2), "sged")
garch21_nor_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,1), "norm")
garch21_std_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,1), "std")
garch21_sstd_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,1), "sstd")
garch21_ged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,1), "ged")
garch21_sged_NI2251 <- fit_gjrgarch(truoc$NI225, arma1NI225, c(2,1), "sged")
#Giai đoạn covid
garch11_nor_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,1), "norm")
garch11_std_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,1), "std")
garch11_sstd_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,1), "sstd")
garch11_ged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,1), "ged")
garch11_sged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,1), "sged")
garch22_nor_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,2), "norm")
garch22_std_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,2), "std")
garch22_sstd_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,2), "sstd")
garch22_ged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,2), "ged")
garch22_sged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,2), "sged")
garch12_nor_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,2), "norm")
garch12_std_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,2), "std")
garch12_sstd_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,2), "sstd")
garch12_ged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,2), "ged")
garch12_sged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(1,2), "sged")
garch21_nor_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,1), "norm")
garch21_std_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,1), "std")
garch21_sstd_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,1), "sstd")
garch21_ged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,1), "ged")
garch21_sged_NI2252 <- fit_gjrgarch(covid$NI225, arma2NI225, c(2,1), "sged")
#Giai đoạn xung đột
garch11_nor_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,1), "norm")
garch11_std_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,1), "std")
garch11_sstd_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,1), "sstd")
garch11_ged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,1), "ged")
garch11_sged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,1), "sged")
garch22_nor_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,2), "norm")
garch22_std_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,2), "std")
garch22_sstd_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,2), "sstd")
garch22_ged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,2), "ged")
garch22_sged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,2), "sged")
garch12_nor_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,2), "norm")
garch12_std_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,2), "std")
garch12_sstd_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,2), "sstd")
garch12_ged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,2), "ged")
garch12_sged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(1,2), "sged")
garch21_nor_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,1), "norm")
garch21_std_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,1), "std")
garch21_sstd_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,1), "sstd")
garch21_ged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,1), "ged")
garch21_sged_NI2253 <- fit_gjrgarch(xung_dot$NI225, arma3NI225, c(2,1), "sged")
#Giai đoạn sau
garch11_nor_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,1), "norm")
garch11_std_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,1), "std")
garch11_sstd_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,1), "sstd")
garch11_ged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,1), "ged")
garch11_sged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,1), "sged")
garch22_nor_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,2), "norm")
garch22_std_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,2), "std")
garch22_sstd_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,2), "sstd")
garch22_ged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,2), "ged")
garch22_sged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,2), "sged")
garch12_nor_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,2), "norm")
garch12_std_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,2), "std")
garch12_sstd_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,2), "sstd")
garch12_ged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,2), "ged")
garch12_sged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(1,2), "sged")
garch21_nor_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,1), "norm")
garch21_std_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,1), "std")
garch21_sstd_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,1), "sstd")
garch21_ged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,1), "ged")
garch21_sged_NI2254 <- fit_gjrgarch(sau$NI225, arma4NI225, c(2,1), "sged")


#Chọn mô hình tối ưu
vni1.model.list <- list(garch11n = garch11_nor_vni1, garch11t = garch11_std_vni1, garch11st = garch11_sstd_vni1, 
                       garch11g = garch11_ged_vni1, garch11sg = garch11_sged_vni1, 
                       garch12n = garch12_nor_vni1, garch12t = garch12_std_vni1, garch12st = garch12_sstd_vni1, 
                       garch12g = garch12_ged_vni1, garch12sg = garch12_sged_vni1,
                       garch21n = garch21_nor_vni1, garch21t = garch21_std_vni1, garch21st = garch21_sstd_vni1, 
                       garch21g = garch21_ged_vni1, garch21sg = garch21_sged_vni1,
                       garch22n = garch22_nor_vni1, garch22t = garch22_std_vni1, garch22st = garch22_sstd_vni1, 
                       garch22g = garch22_ged_vni1, garch22sg = garch22_sged_vni1
                       )

vni1.info.mat <- sapply(vni1.model.list, infocriteria)

rownames(vni1.info.mat) <- rownames(infocriteria(garch11_nor_vni1))

vni1.inds <- which(vni1.info.mat == min(vni1.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.vni1 <- colnames(vni1.info.mat)[vni1.inds[,2]]

vni2.model.list <- list(garch11n = garch11_nor_vni2, garch11t = garch11_std_vni2, garch11st = garch11_sstd_vni2, 
                      garch11sg = garch11_sged_vni2, 
                       garch12n = garch12_nor_vni2, garch12t = garch12_std_vni2, garch12st = garch12_sstd_vni2, 
                       garch12g = garch12_ged_vni2, garch12sg = garch12_sged_vni2,
                       garch21n = garch21_nor_vni2, garch21t = garch21_std_vni2, garch21st = garch21_sstd_vni2, 
                       garch21g = garch21_ged_vni2, garch21sg = garch21_sged_vni2,
                       garch22n = garch22_nor_vni2, garch22t = garch22_std_vni2, garch22st = garch22_sstd_vni2, 
                       garch22g = garch22_ged_vni2, garch22sg = garch22_sged_vni2
                       )

vni2.info.mat <- sapply(vni2.model.list, infocriteria)

rownames(vni2.info.mat) <- rownames(infocriteria(garch11_nor_vni2))

vni2.inds <- which(vni2.info.mat == min(vni2.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.vni2 <- colnames(vni2.info.mat)[vni2.inds[,2]]

vni3.model.list <- list(garch11n = garch11_nor_vni3,  garch11st = garch11_sstd_vni3, 
                       garch11g = garch11_ged_vni3, garch11sg = garch11_sged_vni3, 
                       garch12n = garch12_nor_vni3, garch12t = garch12_std_vni3, garch12st = garch12_sstd_vni3, 
                       garch12g = garch12_ged_vni3, garch12sg = garch12_sged_vni3,
                       garch21n = garch21_nor_vni3, garch21t = garch21_std_vni3, garch21st = garch21_sstd_vni3, 
                       garch21g = garch21_ged_vni3, garch21sg = garch21_sged_vni3,
                       garch22n = garch22_nor_vni3, garch22t = garch22_std_vni3, garch22st = garch22_sstd_vni3, 
                       garch22g = garch22_ged_vni3, garch22sg = garch22_sged_vni3
                       )

vni3.info.mat <- sapply(vni3.model.list, infocriteria)

rownames(vni3.info.mat) <- rownames(infocriteria(garch11_nor_vni3))

vni3.inds <- which(vni3.info.mat == min(vni3.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.vni3 <- colnames(vni3.info.mat)[vni3.inds[,2]]

vni4.model.list <- list(garch11n = garch11_nor_vni4, garch11t = garch11_std_vni4, garch11st = garch11_sstd_vni4, 
                       garch11g = garch11_ged_vni4, garch11sg = garch11_sged_vni4, 
                       garch12n = garch12_nor_vni4, garch12t = garch12_std_vni4, garch12st = garch12_sstd_vni4, 
                       garch12g = garch12_ged_vni4, garch12sg = garch12_sged_vni4
                       )

vni4.info.mat <- sapply(vni4.model.list, infocriteria)

rownames(vni4.info.mat) <- rownames(infocriteria(garch11_nor_vni4))

vni4.inds <- which(vni4.info.mat == min(vni4.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.vni4 <- colnames(vni4.info.mat)[vni4.inds[,2]]

#Cho NI225

NI2251.model.list <- list(garch11n = garch11_nor_NI2251, garch11t = garch11_std_NI2251, garch11st = garch11_sstd_NI2251, 
                       garch11g = garch11_ged_NI2251, garch11sg = garch11_sged_NI2251, 
                       garch12n = garch12_nor_NI2251, garch12t = garch12_std_NI2251, garch12st = garch12_sstd_NI2251, 
                       garch12g = garch12_ged_NI2251, garch12sg = garch12_sged_NI2251,
                       garch21n = garch21_nor_NI2251, garch21t = garch21_std_NI2251, garch21st = garch21_sstd_NI2251, 
                       garch21g = garch21_ged_NI2251, garch21sg = garch21_sged_NI2251,
                       garch22n = garch22_nor_NI2251, garch22t = garch22_std_NI2251, garch22st = garch22_sstd_NI2251, 
                       garch22g = garch22_ged_NI2251, garch22sg = garch22_sged_NI2251
                       )

NI2251.info.mat <- sapply(NI2251.model.list, infocriteria)

rownames(NI2251.info.mat) <- rownames(infocriteria(garch11_nor_NI2251))

NI2251.inds <- which(NI2251.info.mat == min(NI2251.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.NI2251 <- colnames(NI2251.info.mat)[NI2251.inds[,2]]

NI2252.model.list <- list(garch11n = garch11_nor_NI2252, garch11t = garch11_std_NI2252, garch11st = garch11_sstd_NI2252, 
                       garch11g = garch11_ged_NI2252, garch11sg = garch11_sged_NI2252, 
                       garch12n = garch12_nor_NI2252, garch12t = garch12_std_NI2252, garch12st = garch12_sstd_NI2252, 
                       garch12g = garch12_ged_NI2252, garch12sg = garch12_sged_NI2252,
                       garch21n = garch21_nor_NI2252, garch21t = garch21_std_NI2252, garch21st = garch21_sstd_NI2252, 
                       garch21g = garch21_ged_NI2252, garch21sg = garch21_sged_NI2252,
                       garch22n = garch22_nor_NI2252, garch22t = garch22_std_NI2252, garch22st = garch22_sstd_NI2252, 
                       garch22g = garch22_ged_NI2252, garch22sg = garch22_sged_NI2252
                       )

NI2252.info.mat <- sapply(NI2252.model.list, infocriteria)

rownames(NI2252.info.mat) <- rownames(infocriteria(garch11_nor_NI2252))

NI2252.inds <- which(NI2252.info.mat == min(NI2252.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.NI2252 <- colnames(NI2252.info.mat)[NI2252.inds[,2]]

NI2253.model.list <- list(garch11n = garch11_nor_NI2253, garch11t = garch11_std_NI2253, garch11st = garch11_sstd_NI2253, 
                       garch11g = garch11_ged_NI2253, garch11sg = garch11_sged_NI2253, 
                       garch12n = garch12_nor_NI2253, garch12t = garch12_std_NI2253, garch12st = garch12_sstd_NI2253, 
                       garch12g = garch12_ged_NI2253, garch12sg = garch12_sged_NI2253,
                       garch21n = garch21_nor_NI2253, garch21t = garch21_std_NI2253, garch21st = garch21_sstd_NI2253, 
                       garch21g = garch21_ged_NI2253, garch21sg = garch21_sged_NI2253
                       )

NI2253.info.mat <- sapply(NI2253.model.list, infocriteria)

rownames(NI2253.info.mat) <- rownames(infocriteria(garch11_nor_NI2253))

NI2253.inds <- which(NI2253.info.mat == min(NI2253.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.NI2253 <- colnames(NI2253.info.mat)[NI2253.inds[,2]]

NI2254.model.list <- list(garch11t = garch11_std_NI2254, garch11st = garch11_sstd_NI2254, 
                       garch11g = garch11_ged_NI2254, garch11sg = garch11_sged_NI2254, 
                       garch21n = garch21_nor_NI2254, garch21t = garch21_std_NI2254, garch21st = garch21_sstd_NI2254, 
                       garch21g = garch21_ged_NI2254, garch21sg = garch21_sged_NI2254,
                       garch22n = garch22_nor_NI2254, garch22t = garch22_std_NI2254, garch22st = garch22_sstd_NI2254, 
                       garch22g = garch22_ged_NI2254, garch22sg = garch22_sged_NI2254
                       )

NI2254.info.mat <- sapply(NI2254.model.list, infocriteria)

rownames(NI2254.info.mat) <- rownames(infocriteria(garch11_std_NI2254))

NI2254.inds <- which(NI2254.info.mat == min(NI2254.info.mat[1,]), arr.ind=TRUE) #Lay chi so AIC nho nhat

model.NI2254 <- colnames(NI2254.info.mat)[NI2254.inds[,2]]

#Trình bày kết quả
mar_model1 <- data.frame(
  rate = c("VNI", "NI225"),
  truoc = c("ARMA(2,0)-GJR-GARCH(1,1)-Skewed Student t", "ARMA(0,0)-GJR-GARCH(1,1)-Student t"),
  covid = c("ARMA(2,0)-GJR-GARCH(2,2)-Skewed Student t", "ARMA(0,2)-GJR-GARCH(1,1)-Student t"),
  xungdot = c("ARMA(2,2)-GJR-GARCH(1,1)-Skewed Ged", "ARMA(2,2)-GJR-GARCH(1,1)-Skewed Student t"),
  sau = c("ARMA(2,0)-GJR-GARCH(2,2)-Skewed Ged", "ARMA(1,0)-GJR-GARCH(1,1)-Student t"))

# Render the table
kable(mar_model1, col.names = c("Chỉ số", "Giai đoạn trước","Giao đoạn Covid","Giai đoạn xung đột","Giai đoạn sau"), 
      caption = "Bảng 10: Mô hình phân phối biên tối ưu cho các giai đoạn", format = "pandoc", 
      table.attr = "style='width:100%;'") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Bảng 10: Mô hình phân phối biên tối ưu cho các giai đoạn

Chỉ số	Giai đoạn trước	Giao đoạn Covid	Giai đoạn xung đột	Giai đoạn sau
VNI	ARMA(2,0)-GJR-GARCH(1,1)-Skewed Student t	ARMA(2,0)-GJR-GARCH(2,2)-Skewed Student t	ARMA(2,2)-GJR-GARCH(1,1)-Skewed Ged	ARMA(2,0)-GJR-GARCH(2,2)-Skewed Ged
NI225	ARMA(0,0)-GJR-GARCH(1,1)-Student t	ARMA(0,2)-GJR-GARCH(1,1)-Student t	ARMA(2,2)-GJR-GARCH(1,1)-Skewed Student t	ARMA(1,0)-GJR-GARCH(1,1)-Student t

Nhìn chung, các chuỗi lợi suất chứng khoán ở cả 4 giai đoạn đều có phân phối đuôi dày hoặc bất đối xứng chứng tỏ rằng phân phối của dữ liệu rất khác so với phân phối chuẩn.

#Giai đoạn trước
vni.res1 <- residuals(garch11_sstd_vni1)/sigma(garch11_sstd_vni1) #standarđ resid

a1 <- fitdist(distribution = "sstd", vni.res1, control = list()) #fit distribution for this resid

u1 <- pdist(distribution = "sstd", q = vni.res1, mu = a1$pars['mu'] , sigma = a1$pars['sigma'],
           skew = a1$pars['skew'], shape = a1$pars['shape']) # make probability integration

# make probability integration for residual of NI225

NI225.res1 <- residuals(garch11_std_NI2251)/sigma(garch11_std_NI2251) #standarđ resid

b1 <- fitdist(distribution = "std", NI225.res1, control = list()) #fit distribution for this resid

v1 <- pdist(distribution = "std", q = NI225.res1, mu = b1$pars['mu'] , sigma = b1$pars['sigma'],
           skew = b1$pars['skew'], shape = b1$pars['shape']) # make probability integration

#Covid
vni.res2 <- residuals(garch11_sstd_vni2)/sigma(garch11_sstd_vni2) #standarđ resid

a2 <- fitdist(distribution = "sstd", vni.res2, control = list()) #fit distribution for this resid

u2 <- pdist(distribution = "sstd", q = vni.res2, mu = a2$pars['mu'] , sigma = a2$pars['sigma'],
           skew = a2$pars['skew'], shape = a2$pars['shape']) # make probability integration

# make probability integration for residual of NI225

NI225.res2 <- residuals(garch11_std_NI2252)/sigma(garch11_std_NI2252) #standarđ resid

b2 <- fitdist(distribution = "std", NI225.res2, control = list()) #fit distribution for this resid

v2 <- pdist(distribution = "std", q = NI225.res2, mu = b2$pars['mu'] , sigma = b2$pars['sigma'],
           skew = b2$pars['skew'], shape = b2$pars['shape']) # make probability integration

#Xung đột
vni.res3 <- residuals(garch11_sged_vni3)/sigma(garch11_sged_vni3) #standarđ resid

a3 <- fitdist(distribution = "sged", vni.res3, control = list()) #fit distribution for this resid

u3 <- pdist(distribution = "sged", q = vni.res3, mu = a3$pars['mu'] , sigma = a3$pars['sigma'],
           skew = a3$pars['skew'], shape = a3$pars['shape']) # make probability integration

# make probability integration for residual of NI225

NI225.res3 <- residuals(garch11_sstd_NI2253)/sigma(garch11_sstd_NI2253) #standarđ resid

b3 <- fitdist(distribution = "sstd", NI225.res3, control = list()) #fit distribution for this resid

v3 <- pdist(distribution = "sstd", q = NI225.res3, mu = b3$pars['mu'] , sigma = b3$pars['sigma'],
           skew = b3$pars['skew'], shape = b3$pars['shape']) # make probability integration

#Sau
vni.res4 <- residuals(garch11_sged_vni4)/sigma(garch11_sged_vni4) #standarđ resid

a4 <- fitdist(distribution = "sged", vni.res4, control = list()) #fit distribution for this resid

u4 <- pdist(distribution = "sged", q = vni.res4, mu = a4$pars['mu'] , sigma = a4$pars['sigma'],
           skew = a4$pars['skew'], shape = a4$pars['shape']) # make probability integration

# make probability integration for residual of NI225

NI225.res4 <- residuals(garch11_std_NI2254)/sigma(garch11_std_NI2254) #standarđ resid

b4 <- fitdist(distribution = "std", NI225.res4, control = list()) #fit distribution for this resid

v4 <- pdist(distribution = "std", q = NI225.res4, mu = b4$pars['mu'] , sigma = b4$pars['sigma'],
           skew = b4$pars['skew'], shape = b4$pars['shape']) # make probability integration

library(VineCopula)
cop1 <- BiCopSelect(u1,v1,selectioncrit = "AIC",method = "mle")
cop2 <- BiCopSelect(u2,v2,selectioncrit = "AIC",method = "mle")
cop3 <- BiCopSelect(u3,v3,selectioncrit = "AIC",method = "mle")
cop4 <- BiCopSelect(u4,v4,selectioncrit = "AIC",method = "mle")
#est
est_cop1 <- BiCopEst(u1,v1,family = cop1$family,method = "mle",max.df = 30)
est_cop2 <- BiCopEst(u2,v2,family = cop2$family,method = "mle",max.df = 30)
est_cop3 <- BiCopEst(u3,v3,family = cop3$family,method = "mle",max.df = 30)
est_cop4 <- BiCopEst(u4,v4,family = cop4$family,method = "mle",max.df = 30)
#Trình bày kết quả

est1 <- data.frame(mqh = c('Giai đoạn trước' ,'Giai đoạn Covid','Xung đột Nga-Ukraine','Giai đoạn sau'),
                    Copula = c('BB7', 'Survival BB7', 'Clayton','Gauss'),
                    Thamso = c(est_cop1$par, est_cop2$par, est_cop3$par, est_cop4$par),
                  Thamso2 = c(est_cop1$par2, est_cop2$par2, est_cop3$par2, est_cop4$par2),
                   duoi = c(est_cop1$tau, est_cop2$tau, est_cop3$tau, est_cop4$tau),
                  duoiduoi = c(est_cop1$taildep$lower, est_cop2$taildep$lower, est_cop3$taildep$lower, est_cop4$taildep$lower),
                  duoitren = c(est_cop1$taildep$upper, est_cop2$taildep$upper, est_cop3$taildep$upper, est_cop4$taildep$upper),
                  aic = c(est_cop1$AIC, est_cop2$AIC, est_cop3$AIC, est_cop4$AIC),
                  bic = c(est_cop1$BIC, est_cop2$BIC, est_cop3$BIC, est_cop4$BIC))

kable(est1, 
  caption = "Bảng 11: Kết quả ước lượng Copula của TSLN VNI và NI225 từng giai đoạn", 
  col.names = c("Chỉ số","Copula", "Par","Par2", "Hệ số Kendall","Đuôi dưới","Đuôi trên","AIC","BIC"), 
  format = 'pandoc') %>% 
  kable_styling(
            bootstrap_options = c("striped", "hover", "condensed"), 
            full_width = F)

Bảng 11: Kết quả ước lượng Copula của TSLN VNI và NI225 từng giai đoạn

Chỉ số	Copula	Par	Par2	Hệ số Kendall	Đuôi dưới	Đuôi trên	AIC	BIC
Giai đoạn trước	BB7	0.3669004	6.249583	0.2391621	0.1080972	0.1080972	-55.85084	-48.00233
Giai đoạn Covid	Survival BB7	0.2708931	11.736857	0.1746380	0.0183581	0.0183581	-33.24308	-25.07864
Xung đột Nga-Ukraine	Clayton	1.2721648	0.000000	0.1340406	0.2756392	0.0000000	-15.35277	-11.99619
Giai đoạn sau	Gauss	1.2055585	0.000000	0.1705089	0.0000000	0.2229416	-31.20045	-27.49334

Bảng trên trình bày kết quả ước lượng các mô hình Copula cho từng giai đoạn khác nhau dựa trên hai chỉ số tỷ suất lợi nhuận của VNI và NI225. Trong giai đoạn trước các cú sốc kinh tế, mô hình BB7 là Copula phù hợp nhất đã cho thấy mối tương quan dương giữa VNI và NI225 với hệ số Kendall là 0.2391621, xác suất xảy ra sự kiện cực đoan ở phần trên (0.1080972) bằng phần dưới (0.1080972) chứng tỏ rằng mức độ phản ứng của TTCK Việt Nam sẽ tương đương khi TTCK Nhật Bản tăng giá cực biên so với việc giảm giá. Khi bước vào giai đoạn Covid, mô hình Survival BB7 phản ánh sự thay đổi trong cấu trúc phụ thuộc, khi cho thấy xác suất xảy ra sự kiện cực đoan ở phần dưới (0.0183581) tương đương trên (0.0183581), biểu hiện tác động mạnh mẽ của đại dịch đã ảnh hưởng đến cả hai TTCK khiến cả hai cùng giảm, Tuy nhiên, hệ số Kendall của giai đoạn này (0.1746380) đã giảm nhẹ so với giai đoạn trước đó. Tiếp đến, trong giai đoạn xung đột Nga-Ukraine, mô hình Clayton được sử dụng để phản ánh sự phụ thuộc cao hơn giữa hai chỉ số, với xác suất xảy ra sự kiện cực đoan chủ yếu tập trung ở phần dưới (0.2756392), biểu hiện tình hình bất ổn do xung đột đã đẩy mức độ nhạy cảm đối với các tin tức tiêu cực của các TTCK lên, điều này chỉ ra rằng nếu một trong hai TTCK sụp đổ thì xác suất TTCK còn lại sụp đổ theo là khoảng 27,5% nhưng nếu một trong hai thị trường tăng giá bùng nổ thì không ảnh hưởng gì đến TTCK còn lại, thêm nữa hệ số Kendall của giai đoạn này (0.1340406) là cao nhất trong các giai đoạn chứng tỏ mức độ nhạy cảm của hai TTCK đã gia tăng trong giai đoạn xảy ra xung đột Nga-Ukraine. Cuối cùng, trong giai đoạn sau cú sốc kinh tế, Copula được sử dụng để mô hình hóa sự phụ thuộc giữa hai thị trường là mô hình Gauss chứng tỏ rằng giai đoạn nghiên cứu này là giai đoạn ổn định và có rất ít các biến động cực biên, hệ số Kendall của giai đoạn này (0.1705089) là thấp nhất trong tất cả các giai đoạn, chỉ ra rằng mối quan hệ giữa hai chỉ số trở nên rất yếu sau những biến động lớn trước đó.

library(VC2copula)
set.seed(123)
#Mô phỏng copula cho giai đoạn trước
library(copula)
cop_1 <- normalCopula(est_cop1$par, dim = 2)

#Mô phỏng copula cho giai đoạn covid
cop_2 <- normalCopula(est_cop2$par, dim = 2)

#Mô phỏng copula cho giai đoạn xung đột
cop_3 <- claytonCopula(est_cop3$par, dim = 2)

#Mô phỏng Copula cho giai đoạn sau
cop_4 <- claytonCopula(est_cop4$par, dim = 2)

per1 <- persp_plot(cop_1, 'cop1.png', '#FFFF00')
per2 <- persp_plot(cop_2, 'cop2.png', '#FF9900')
per3 <- persp_plot(cop_3, 'cop3.png', '#000080')
per4 <- persp_plot(cop_4, 'cop4.png', '#FF3300')


legend <- legendGrob(
  labels = c("Trước - BB7", "Covid - Sur BB7","Xung đột - Clayton", "Sau - Gauss"), pch = 15,
  gp = gpar(col = c('#FFFF00', '#FF9900', '#000080','#FF3300'), fill = c('#FFFF00', '#FF9900', '#000080','#FF3300'))
)

grid.arrange( per1, per2, per3,per4, legend, ncol = 3, 
  layout_matrix = rbind(c(1, 2, 5), c(3,4,5)),
  widths = c(2, 2, 1), 
  top = textGrob("Hình 14: Phối cảnh PDF của Copula các giai đoạn", 
                 gp = gpar(fontsize = 15, font = 2))
             )

CHƯƠNG 5: KẾT LUẬN VÀ KHUYẾN NGHỊ

5.1 Kết luận

Bài tiểu luận đã tiến hành nghiên cứu phản ứng của thị trường chứng khoán (TTCK) Việt Nam trước sự biến động của TTCK Nhật Bản, thông qua việc phân tích chuỗi tỷ suất sinh lợi của hai chỉ số tiêu biểu là VNINDEX và NI225. Cụ thể, nghiên cứu đã xem xét cấu trúc phụ thuộc giữa hai TTCK trong giai đoạn từ đầu năm 2018 đến đầu tháng 7/2024 bằng phương pháp Copula có điều kiện. Nghiên cứu cũng đã áp dụng phương pháp này cho từng giai đoạn riêng biệt để đánh giá mức độ phản ứng của TTCK Việt Nam qua các thời kỳ ổn định và biến động mạnh. Phương pháp copula là một công cụ mạnh mẽ trong việc phân tích cấu trúc phụ thuộc giữa các biến ngẫu nhiên thông qua xây dựng hàm phân phối đa biến, đồng thời cho phép mô hình hóa sự phụ thuộc đuôi giữa hai chuỗi tỷ suất sinh lợi trên thị trường tài chính, đặc biệt khi các phân phối đơn biến phức tạp. Một số mô hình phân tích thông thường giả định rằng chuỗi lợi suất tuân theo phân phối chuẩn, tuy nhiên, thực tế cho thấy các chuỗi này có sự biến động lớn và ngẫu nhiên. Do đó, phương pháp copula, thông qua các phân phối biên khác nhau, cho phép phân tích sự phụ thuộc đối xứng và bất đối xứng. Nghiên cứu này đã áp dụng nhiều hàm copula như Gauss, Student-t, Clayton, Rotated-Clayton, BB1, BB7 để đánh giá mức độ phụ thuộc giữa TTCK Việt Nam và TTCK Nhật Bản.

Kết quả mô hình Copula đã chứng tỏ được sự tồn tại của mối tương quan tích cực giữa TTCK Việt Nam đối với TTCK Nhật Bản trong giai đoạn đầu năm 2018 đến tháng 7/2024. Mô hình Copula phù hợp nhất cho dữ liệu là BB7với hệ số Kendall là 0.2391621 cho thấy sự phụ thuộc không đáng kể. Bên cạnh đó, kết quả còn cho thấy sự phụ thuộc đuôi dưới (0.0.1080972) là nhỏ hơn phụ thuộc đuôi trên (0.1080972), điều này chỉ ra rằng xác suất để hai TTCK bùng nổ cùng nhau là lớn hơn xác suất cả hai TTCK cùng sụp đổ. Ngoài ra, nghiên cứu cũng đã chứng tỏ được mối tương quan của hai TTCK gia tăng trong thời kỳ xảy ra cú sốc kinh tế (Xung đột Nga-Ukraine) với hệ số phụ thuộc lớn nhất trong 4 giai đoạn nghiên cứu. Điều này lần nữa khẳng định rằng, khi xảy ra các cú sốc kinh tế, mối tương quan của hai TTCK sẽ gia tăng và khả năng hai TTCK cùng giảm giá cực biên cũng gia tăng đáng kể.

5.2 Khuyến nghị

Dựa trên kết quả từ mô hình copula BB7 về mối quan hệ giữa tỷ suất lợi nhuận của chỉ số VNI và NI225 từ đầu năm 2018 đến tháng 7/2024, có thể rút ra một số khuyến nghị quan trọng cho các nhà đầu tư và nhà hoạch định chính sách. Đối với các nhà đầu tư, việc đa dạng hóa danh mục đầu tư là rất cần thiết do mức độ tương quan dương nhưng không quá mạnh giữa VNI và NI225, giúp giảm thiểu rủi ro tổng thể. Kết quả cũng cho thấy sự phụ thuộc gia tăng trong các giai đoạn biến động mạnh, đặc biệt là trong các giai đoạn có sự biến động cao như đại dịch Covid-19 và xung đột Nga-Ukraine. Do đó, nhà đầu tư cần thường xuyên cập nhật thông tin về kinh tế và chính trị để tránh những rủi ro bất lợi. Việc theo dõi sát sao các sự kiện này giúp nhà đầu tư có thể điều chỉnh chiến lược đầu tư kịp thời. Nhà đầu tư cũng cần kết hợp phân tích kỹ thuật và phân tích cơ bản, cùng với việc sử dụng các công cụ phái sinh để bảo vệ danh mục đầu tư khỏi các sự kiện bất ngờ. Ngoài ra, quản lý rủi ro một cách cẩn thận thông qua các công cụ phòng ngừa rủi ro là rất quan trọng, bởi có khả năng xảy ra các biến động cực đoan cùng chiều, đặc biệt khi thị trường tăng mạnh. Nhà đầu tư cũng nên tận dụng cơ hội khi một trong hai chỉ số có dấu hiệu tăng trưởng mạnh, nhờ vào mức độ phụ thuộc đuôi trên cao hơn.

Đối với các nhà hoạch định chính sách, việc giám sát chặt chẽ các yếu tố kinh tế và chính trị ảnh hưởng đến cả hai thị trường là rất quan trọng. Cần duy trì các chính sách kinh tế ổn định để tạo ra một môi trường đầu tư an toàn và giảm thiểu các biến động cực đoan. Đồng thời, cần xây dựng các biện pháp hỗ trợ kinh tế trong thời kỳ khủng hoảng nhằm tăng cường sự ổn định kinh tế. Các chính sách quản lý rủi ro tài chính cũng cần được thực hiện một cách chặt chẽ để ngăn chặn sự lan truyền rủi ro hệ thống. Hơn nữa, việc hợp tác quốc tế giữa các cơ quan quản lý tài chính có thể giúp đồng bộ hóa các biện pháp phòng ngừa rủi ro và đảm bảo sự ổn định của thị trường. Các nhà hoạch định chính sách cũng cần đảm bảo tính minh bạch và cung cấp thông tin kịp thời, giúp nhà đầu tư đưa ra các quyết định chính xác hơn, từ đó ổn định thị trường. Ngoài ra, khuyến khích các khoản đầu tư bền vững và dài hạn sẽ giúp giảm thiểu rủi ro và hỗ trợ sự phát triển bền vững của nền kinh tế trong dài hạn.

TÀI LIỆU THAM KHẢO

1. Yasushi Hamao, Ronald W. Masulis, Victor Ng, Correlations in Price Changes and Volatility across International Stock Markets, The Review of Financial Studies, Volume 3, Issue 2, April 1990, Pages 281–307

2. Longin, F., & Solnik, B. (1995). Is the correlation in international equity returns constant: 1960–1990?. Journal of international money and finance, 14(1), 3-26.

3. Xu, H., & Hamori, S. (2012). Dynamic linkages of stock prices between the BRICs and the United States: Effects of the 2008–09 financial crisis. Journal of Asian Economics, 23(4), 344-352.

4. Gupta, R., & Guidi, F. (2012). Cointegration relationship and time varying co movements among Indian and Asian developed stock markets. International Review of Financial Analysis, 21, 10-22.

5. Karanasos, M., Yfanti, S., & Karoglou, M. (2016). Multivariate FIAPARCH modelling of financial markets with dynamic correlations in times of crisis. International Review of Financial Analysis, 45, 332-349.

6. Poon, S. H., Rockinger, M., & Tawn, J. (2004). Extreme value dependence in financial markets: Diagnostics, models, and financial implications. The Review of Financial Studies, 17(2), 581-610.

7. Basher, S. A., Nechi, S., & Zhu, H. (2014). Dependence patterns across Gulf Arab stock markets: A copula approach. Journal of Multinational Financial Management, 25, 30-50.

8. Jondeau, E., & Rockinger, M. (2006). The copula-garch model of conditional dependencies: An international stock market application. Journal of international money and finance, 25(5), 827-853.

9. Mensah, J. O., & Alagidede, P. (2017). How are Africa’s emerging stock markets related to advanced markets? Evidence from copulas. Economic Modelling, 60, 1-10.

10. Hussain, S. I., & Li, S. (2018). The dependence structure between Chinese and other major stock markets using extreme values and copulas. International Review of Economics & Finance, 56, 421-437.

11. Le Van Thu (2022), PHẢN ỨNG CỦA THỊ TRƯỜNG CHỨNG KHOÁN VIỆT NAM ĐỐI VỚI BIẾN ĐỘNG CỦA THỊ TRƯỜNG CHỨNG KHOÁN QUỐC TẾ . Bollerslev, T. (1990). Aconditionally Heteroskedastic Time Series Model for Speculative. Tim.

12. Bollerslev, T. (1990). Modelling The Coherence in Short-run Nominal Exchange Rates A Multivariate Generalized ARCH Model. Tim.

13. Hien, K. a. (2003). A Study on Technical Efficiency of Rice Production in the Mekong Delta, Vietnam by Stochastic Frontier Analysis”. Faculty of Agriculture, Kyushu University.

14. Li, W. K. (1994). On The Squared Residual Autocorrelations in Non-linear Time series with Conditional Heteroskedasticity. Wai Keung Li.

15. Mika Meitz, P. S. (2008). Parameter Estimation in Nonlinear AR-GARCH Models. Mika Meitz.

16. Phuc . (2022). “Quantifying Heterogeneity, Heteroscedasticity, and Publication Bias Effects on Technical Efficiency Estimates of Rice Farming:A Meta‐Regression Analysis”. Journal of Agricultural Economics.

17. Sadorsky, P. (2006). Modeling and forestcasting petroleum futures volatility. In P. Sadorsky, Energy economics (pp. 467 - 468).