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
| 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 = 'DataMidas.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_H = DATA_H(1:end,2:r+1); % variables explicatives mensuelles
%% 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]); |
