clc;clear;


% load climate and economic data
load('climate_data_obs.mat'); % nino3.4 index for 1950-2021
load('gdp.mat'); % gdp data for 1960-2019

% fit econometric model for all countries

delta_gdp_allcountry=diff(log(gdp_pc_country(:,1:60)),1,2);
% calculate linear and quadratic time trend in each model
for i=1:size(delta_gdp_allcountry,1)
    t=1961:1:2019;
    t(isnan(delta_gdp_allcountry(i,:)))=nan;
    p=polyfit(rmmissing(t),rmmissing(delta_gdp_allcountry(i,:)),2);
    theta_qua_allcountry(i)=p(1);
    theta_allcountry(i)=p(2);
end
[djf_nino34_hadisst_m,~]=meshgrid(djf_nino34_hadisst_nor(1:end-1),psi_country_nino34);
nino34_allcountry=djf_nino34_hadisst_m;
nino34qua_allcountry=djf_nino34_hadisst_m.^2;
% test lag 5
nino34_0_allcountry=nino34_allcountry(:,13:end);
nino34_1_allcountry=nino34_allcountry(:,12:end-1);
nino34_2_allcountry=nino34_allcountry(:,11:end-2);
nino34_3_allcountry=nino34_allcountry(:,10:end-3);
nino34_4_allcountry=nino34_allcountry(:,9:end-4);
% quadratic term
nino34qua_0_allcountry=nino34qua_allcountry(:,13:end);
nino34qua_1_allcountry=nino34qua_allcountry(:,12:end-1);
nino34qua_2_allcountry=nino34qua_allcountry(:,11:end-2);
nino34qua_3_allcountry=nino34qua_allcountry(:,10:end-3);
nino34qua_4_allcountry=nino34qua_allcountry(:,9:end-4);
nino34qua_5_allcountry=nino34qua_allcountry(:,8:end-5);
% mean temperature and precipitation
tas0_allcountry=annual_tas_era5_popweight_country(:,12:end-2);
tas1_allcountry=annual_tas_era5_popweight_country(:,11:end-3);
tas2_allcountry=annual_tas_era5_popweight_country(:,10:end-4);
tas3_allcountry=annual_tas_era5_popweight_country(:,9:end-5);
tas4_allcountry=annual_tas_era5_popweight_country(:,8:end-6);
tas5_allcountry=annual_tas_era5_popweight_country(:,7:end-7);
pr0_allcountry=annual_pr_era5_popweight_country(:,12:end-2);
pr1_allcountry=annual_pr_era5_popweight_country(:,11:end-3);
pr2_allcountry=annual_pr_era5_popweight_country(:,10:end-4);
pr3_allcountry=annual_pr_era5_popweight_country(:,9:end-5);
pr4_allcountry=annual_pr_era5_popweight_country(:,8:end-6);
pr5_allcountry=annual_pr_era5_popweight_country(:,7:end-7);
year_used=1961:2019;

delta_gdp_allcountry=delta_gdp_allcountry(:);% gdp change
nino34_0_allcountry=nino34_0_allcountry(:); 
nino34_1_allcountry=nino34_1_allcountry(:); 
nino34_2_allcountry=nino34_2_allcountry(:); 
nino34_3_allcountry=nino34_3_allcountry(:); 
nino34_4_allcountry=nino34_4_allcountry(:); 
nino34_5_allcountry=nino34_5_allcountry(:); 
nino34qua_0_allcountry=nino34qua_0_allcountry(:); 
nino34qua_1_allcountry=nino34qua_1_allcountry(:); 
nino34qua_2_allcountry=nino34qua_2_allcountry(:); 
nino34qua_3_allcountry=nino34qua_3_allcountry(:); 
nino34qua_4_allcountry=nino34qua_4_allcountry(:); 
nino34qua_5_allcountry=nino34qua_5_allcountry(:); 
tas0_allcountry=tas0_allcountry(:);
tas1_allcountry=tas1_allcountry(:);
tas2_allcountry=tas2_allcountry(:);
tas3_allcountry=tas3_allcountry(:);
tas4_allcountry=tas4_allcountry(:);
tas5_allcountry=tas5_allcountry(:);
tas_qua0_allcountry=tas0_allcountry.^2;
tas_qua1_allcountry=tas1_allcountry.^2;
tas_qua2_allcountry=tas2_allcountry.^2;
tas_qua3_allcountry=tas3_allcountry.^2;
tas_qua4_allcountry=tas4_allcountry.^2;
tas_qua5_allcountry=tas5_allcountry.^2;
pr0_allcountry=pr0_allcountry(:);
pr1_allcountry=pr1_allcountry(:);
pr2_allcountry=pr2_allcountry(:);
pr3_allcountry=pr3_allcountry(:);
pr4_allcountry=pr4_allcountry(:);
pr5_allcountry=pr5_allcountry(:);
pr_qua0_allcountry=pr0_allcountry.^2;
pr_qua1_allcountry=pr1_allcountry.^2;
pr_qua2_allcountry=pr2_allcountry.^2;
pr_qua3_allcountry=pr3_allcountry.^2;
pr_qua4_allcountry=pr4_allcountry.^2;
pr_qua5_allcountry=pr5_allcountry.^2;
country=categorical(repmat(country_list,59,1));% country fixed effect
theta_allcountry=repmat(theta_allcountry',59,1); % linear time trend
theta_qua_allcountry=repmat(theta_qua_allcountry',59,1); % quatratic time trend
year_used=sort(repmat(year_used(:),length(country_list),1));
year_qua=year_used.^2;
trend_allcountry=theta_allcountry.*year_used;
trend_qua_allcountry=theta_qua_allcountry.*year_qua;

% combine all data incorporated in economtreic model
table_data_nino34_allcountry=table(delta_gdp_allcountry,delta_gdp_allcountry_agri,delta_gdp_allcountry_nonagri,...
    nino34_0_allcountry,nino34_1_allcountry,nino34_2_allcountry,nino34_3_allcountry,nino34_4_allcountry,nino34_5_allcountry,...
    nino34qua_0_allcountry,nino34qua_1_allcountry,nino34qua_2_allcountry,nino34qua_3_allcountry,nino34qua_4_allcountry,nino34qua_5_allcountry,...
    tas0_allcountry,tas1_allcountry,tas2_allcountry,tas3_allcountry,tas4_allcountry,tas5_allcountry,...
    tas_qua0_allcountry,tas_qua1_allcountry,tas_qua2_allcountry,tas_qua3_allcountry,tas_qua4_allcountry,tas_qua5_allcountry,...
    pr0_allcountry,pr1_allcountry,pr2_allcountry,pr3_allcountry,pr4_allcountry,pr5_allcountry,...
    pr_qua0_allcountry,pr_qua1_allcountry,pr_qua2_allcountry,pr_qua3_allcountry,pr_qua4_allcountry,...
    country,trend_allcountry,trend_qua_allcountry);

% fit model
mdl_nino34_lag0_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34qua_0_allcountry+tas0_allcountry+tas_qua0_allcountry+pr0_allcountry+pr_qua0_allcountry+country+trend_allcountry+trend_qua_allcountry');
mdl_nino34_lag1_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+tas0_allcountry+tas1_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+pr0_allcountry+pr1_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+country+trend_allcountry+trend_qua_allcountry');
mdl_nino34_lag2_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+country+trend_allcountry+trend_qua_allcountry');
mdl_nino34_lag3_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+country+trend_allcountry+trend_qua_allcountry');
mdl_nino34_lag4_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34_4_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+nino34qua_4_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas4_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+tas_qua4_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr4_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+pr_qua4_allcountry+country+trend_allcountry+trend_qua_allcountry');
mdl_nino34_lag5_allcountry=fitlm(table_data_nino34_allcountry,...
    'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34_4_allcountry+nino34_5_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+nino34qua_4_allcountry+nino34qua_5_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas4_allcountry+tas5_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+tas_qua4_allcountry+tas_qua5_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr4_allcountry+pr5_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+pr_qua4_allcountry+pr_qua5_allcountry+country+trend_allcountry+trend_qua_allcountry');

% extract coefficient of each term
coef_nino34_lag0_allcountry=mdl_nino34_lag0_allcountry.Coefficients.Estimate(1:end);
coef_nino34_lag1_allcountry=mdl_nino34_lag1_allcountry.Coefficients.Estimate(1:end);
coef_nino34_lag2_allcountry=mdl_nino34_lag2_allcountry.Coefficients.Estimate(1:end);
coef_nino34_lag3_allcountry=mdl_nino34_lag3_allcountry.Coefficients.Estimate(1:end);
coef_nino34_lag4_allcountry=mdl_nino34_lag4_allcountry.Coefficients.Estimate(1:end);
coef_nino34_lag5_allcountry=mdl_nino34_lag5_allcountry.Coefficients.Estimate(1:end);

% bootstrap
rng(123);
for i=1:1000
    country_here=randi(length(country_list),[1 length(country_list)]);
    year_here=1961:2019;
    comb_here=[repmat(country_here(:),59,1) sort(repmat(year_here(:),length(country_list),1))];
    [~,ind]=ismember(comb_here,[repmat((1:length(country_list))',59,1) table_data_nino34_allcountry.year_used],'rows');
    table_data_here=table_data_nino34_allcountry(ind,:);

    mdl_nino34_lag0_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34qua_0_allcountry+tas0_allcountry+tas_qua0_allcountry+pr0_allcountry+pr_qua0_allcountry+country+trend_allcountry+trend_qua_allcountry');
    mdl_nino34_lag1_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+tas0_allcountry+tas1_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+pr0_allcountry+pr1_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+country+trend_allcountry+trend_qua_allcountry');
    mdl_nino34_lag2_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+country+trend_allcountry+trend_qua_allcountry');
    mdl_nino34_lag3_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+country+trend_allcountry+trend_qua_allcountry');
    mdl_nino34_lag4_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34_4_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+nino34qua_4_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas4_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+tas_qua4_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr4_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+pr_qua4_allcountry+country+trend_allcountry+trend_qua_allcountry');
    mdl_nino34_lag5_allcountry_bs_here=fitlm(table_data_here,...
        'delta_gdp_allcountry~nino34_0_allcountry+nino34_1_allcountry+nino34_2_allcountry+nino34_3_allcountry+nino34_4_allcountry+nino34_5_allcountry+nino34qua_0_allcountry+nino34qua_1_allcountry+nino34qua_2_allcountry+nino34qua_3_allcountry+nino34qua_4_allcountry+nino34qua_5_allcountry+tas0_allcountry+tas1_allcountry+tas2_allcountry+tas3_allcountry+tas4_allcountry+tas5_allcountry+tas_qua0_allcountry+tas_qua1_allcountry+tas_qua2_allcountry+tas_qua3_allcountry+tas_qua4_allcountry+tas_qua5_allcountry+pr0_allcountry+pr1_allcountry+pr2_allcountry+pr3_allcountry+pr4_allcountry+pr5_allcountry+pr_qua0_allcountry+pr_qua1_allcountry+pr_qua2_allcountry+pr_qua3_allcountry+pr_qua4_allcountry+pr_qua5_allcountry+country+trend_allcountry+trend_qua_allcountry');
    coef_nino34_lag0_allcountry_bs(:,i)=mdl_nino34_lag0_allcountry_bs_here.Coefficients.Estimate(1:end);
    coef_nino34_lag1_allcountry_bs(:,i)=mdl_nino34_lag1_allcountry_bs_here.Coefficients.Estimate(1:end);
    coef_nino34_lag2_allcountry_bs(:,i)=mdl_nino34_lag2_allcountry_bs_here.Coefficients.Estimate(1:end);
    coef_nino34_lag3_allcountry_bs(:,i)=mdl_nino34_lag3_allcountry_bs_here.Coefficients.Estimate(1:end);
    coef_nino34_lag4_allcountry_bs(:,i)=mdl_nino34_lag4_allcountry_bs_here.Coefficients.Estimate(1:end);
    coef_nino34_lag5_allcountry_bs(:,i)=mdl_nino34_lag5_allcountry_bs_here.Coefficients.Estimate(1:end);
end