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]); |
Partager