1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
| %--------------------------------------------------------------------------
%% OBJET: Estimation d'un modèle MIDAS avec une fonction de poids almon exp
% Application à la croissance trimestrielle du PIB
%--------------------------------------------------------------------------
% DESCRIPTIF:
% Programme principal:
% 1) charge les données de haute et basse fréquences
% 2) appelle l'algorithme d'estimation par MCNL
% 3) affiche les résultats d'estimation
% 4) représente les fonctions de poids
% Fonction de poids almonf
%--------------------------------------------------------------------------
clear; clc; close all;
%% Paramètres à modifier
r = 3; % nb de variables explicatives mensuelles
K = 12; % nb de retards du polynome Almon (retard midas)
p_ar = 1; % nb de termes autorégressifs
bounds = 1; % si 1, on impose theta2<0 dans la fonction almon
delai_m = 0; % si >0, on ne connait pas la variable mensuelle sur les delai_m derniers mois du trimestre
options = optimset('TolX',10^(-10),'TolFun',10^(-10),'MaxIter',1000,'MaxFunEvals',10^5,'Display','off');
% Options de l'Optimisation
% conditions initiales
param_init = [1;-1;zeros(p_ar,1);1;0.5;-0.15] + randn(5+p_ar,1);
% paramètres initiaux : c, delta, phi, beta_1 ... beta_r
% theta1_1 ... theta1_r theta2_1... theta2_r
%% Chargement de la BDD
file = 'DATA.xlsx';
DATA_L = xlsread(file,1); % endog = variable expliquée trimestrielle
date_L = DATA_L(1:end,1); % les dates trimestrielles en format numérique
endog_L = DATA_L(1:end,2); % le PIB trimestriel
indic_L = DATA_L(1:end,3); % l'indicatrice du trimestre 2009T1
DATA_H = xlsread(file,2); % variables explicatives mensuelles
date_H = DATA_H(1:end,1); % les dates mensuelles en format numérique
exog = DATA_H(1:end,2:r+1); % variables explicatives mensuelles
%Retraitement des données
oil = diff(log(exog(:,1)));
stock = diff(log(exog(:,2)));
spread = diff(exog(:,3));
exog_H = [oil stock spread];
%% Construction des variables expliquée et explicatives
% Termes autorégressif
X_ar = zeros(size(endog_L,1),p_ar);
for i_p = 1:p_ar
X_ar(:,i_p) = [repmat(mean(endog_L),i_p,1);endog_L(1:end-i_p)];
end
% Recherche des dates communes aux vecteurs H et L
[date_com,index_L,index_H]=intersect(date_L,date_H);
% date_com = dates communes des 2 vecteurs,
% index_L = indice des dates dans le vecteur trimestriel (en principe, 1,2,3...,T)
% index_H = indice des dates dans le vecteur mensuel
% variable expliquée de fréquence trimestrielle
Y = endog_L(index_L);
% variables explicatives de fréquence mensuelle
x_midas_j = zeros(length(index_H),K);
for j=1:length(index_H)
x_midas_j(j,:)=(index_H(j)-delai_m:-1:index_H(j)-K+1-delai_m);
end
% Calcul du nb de valeurs manquantes à retirer
if sum(sum(x_midas_j<=0))>0
[a,b]=find(x_midas_j<=0);
a = unique(a); % lignes à supprimer car comportant indices <0
x_midas_j(x_midas_j<=0)=1;
else
a=0;
end
a = max([a;p_ar]); % termes autorégressifs
X_midas = [];
for i=1:r
temp = exog_H(:,i);
X_midas = [X_midas temp(x_midas_j)];
end
% matrice des explicatives par bloc de K colonnes
% colonnes 1:K : retards de la var 1, K+1:2K : var 2 etc
% on supprime les premiers points du fait des K retards
Y = Y(max(a)+1:end,:);
X_midas = X_midas(max(a)+1:end,:);
X_control = [indic_L(max(a)+1:end,:) X_ar(max(a)+1:end,:)];
n_X_H = size(X_midas,2)/K;
n_X_L = size(X_control,2);
T = size(Y,1);
%% Estimation du modèle
[paramopt,res,code,tstat] = estim_midas(Y,X_midas,X_control,param_init,bounds,options);
%% Tests de diagnostic
[~,pval_JB] = jbtest(res); % JB test
[~,pval_LB] = lbqtest(res,'lags',[4,12,20]); % Ljung Box test
[~,pval_ARCH] = lbqtest(res.^2,'lags',[4,12,20]); % McLeod-Li test
%% Affichage des résultats d'estimation
format compact
disp('Résultats d''estimation');
disp(['Exit condition: ' num2str(code)]);
% >0: Function converged to a solution x
% 0: The maximum number of function evaluations or iterations was exceeded
% <0: The function did not converge to a solution
disp(['R-squared: ' num2str(1-(res'*res)/sum((Y-mean(Y))'*(Y-mean(Y))))]);
disp(['Jarque Bera test: ' num2str(pval_JB)]);
disp(['Ljung Box test: ' num2str(pval_LB)]);
disp(['ARCH effects test: ' num2str(pval_ARCH)]);
nom_coef = ['cste ';'I_09T1';'phi ';'beta ';'theta1';'theta2'];
tab_resu = table(nom_coef,round(paramopt,3),round(tstat,3),'VariableNames',{'Coef','bet','tstat'});
disp(tab_resu)
% Affichage des poids
theta= paramopt(end-1:end)';
poids = poids_almon(theta,K); % Kxr
x = (0:11)';
bar(x,poids); xlim([0 11]); |
