step0_source_data.m

%% Step 0: Fetch source data (OWID + World Bank) and build analysis-ready transforms
% Countries: GRC (Greece), BEL (Belgium)
% Outputs (under raw_data):
%   - OWID_share-electricity-renewables_ALL.csv
%   - WB_<INDICATOR>_<CTY>.csv                (3-col: Year, Value)
%   - WB_<INDICATOR>_<CTY>_META.txt           (provenance)
%   - merged_<CTY>_raw.csv                    (Year, REN_share_pct, GDP_2015USD, K_2015USD, X_2015USD, IMP_2015USD)
%   - merged_<CTY>_META.txt                   (provenance)
%   - merged_<CTY>_transforms.csv             (Year, gdp, k, nx, ren, ren_alt_log)
%   - merged_<CTY>_TRANSFORMS.txt             (units and transform contract)
%
% Transform meanings for VAR in differences:
%   gdp  = Δln(GDP_2015USD)
%   k    = Δln(K_2015USD)                 [GFCF]
%   nx   = Δ(NX/GDP)*100                  [percentage points; NX = X − IMP]
%   ren  = Δlogit(REN_share)              [REN_share = REN_share_pct/100 ∈ (0,1)]
%   ren_alt_log = Δln(REN_share)          [robustness only; requires REN_share>0]
%
% Notes:
%   - OWID “share-electricity-renewables” is fetched via the stable Chart CSV endpoint;
%     we filter GRC/BEL locally.
%   - World Bank indicators are constant 2015 US$:
%       GDP: NY.GDP.MKTP.KD
%       K  : NE.GDI.FTOT.KD
%       X  : NE.EXP.GNFS.KD
%       IMP: NE.IMP.GNFS.KD
%   - Pagination is handled robustly; we try ISO-3 then ISO-2 codes for WB.

clear; clc;

%% ---- Config ----
countries = ["GRC","BEL"];              % ISO-3 (OWID)
analysisWindow = [1990 2023];           % inclusive
baseDir   = "C:\Users\User\OneDrive\Υπολογιστής\Strave\Strave_MATLAB";
rawDir    = fullfile(baseDir, "raw_data");
if ~exist(rawDir, 'dir'); mkdir(rawDir); end

% World Bank indicators (constant 2015 US$)
WB = struct( ...
  'GDP', "NY.GDP.MKTP.KD", ...        % GDP (constant 2015 US$)
  'K',   "NE.GDI.FTOT.KD", ...        % Gross fixed capital formation (constant 2015 US$)
  'X',   "NE.EXP.GNFS.KD", ...        % Exports of goods & services (constant 2015 US$)
  'IMP', "NE.IMP.GNFS.KD" ...         % Imports of goods & services (constant 2015 US$)
);

%% ---- 0.1 Download OWID chart CSV (stable endpoint) ----
% Canonical Chart CSV: columns include Entity, Code, Year, share_electricity_renewables
owidUrl = "https://ourworldindata.org/grapher/share-electricity-renewables.csv";
owidCsvPath = fullfile(rawDir, "OWID_share-electricity-renewables_ALL.csv");
fprintf("Downloading OWID renewables CSV …\n");
try
    websave(owidCsvPath, owidUrl);
    fprintf("Saved: %s\n", owidCsvPath);
catch ME
    error("Failed to fetch OWID CSV: %s", ME.message);
end

% Write minimal provenance
fid = fopen(fullfile(rawDir, "OWID_share-electricity-renewables_META.txt"), 'w');
fprintf(fid, "Source: Our World in Data (Chart CSV)\nURL: %s\nFetched: %s\n", owidUrl, datestr(now,'yyyy-mm-dd HH:MM:SS'));
fprintf(fid, "Series: share_electricity_renewables (percent, 0–100)\nLicense: CC BY (per OWID; underlying sources may vary)\n");
fclose(fid);

% Read OWID with original header names preserved
T_owid_all = readtable(owidCsvPath, 'VariableNamingRule','preserve');
have = T_owid_all.Properties.VariableNames;

% Resolve Code and Year columns
codeIdx = find(strcmpi(have,'Code') | strcmpi(have,'iso_code'), 1);
yearIdx = find(strcmpi(have,'Year'), 1);

% Resolve the renewables share column (try known names; otherwise fallback)
valCandidates = {'share_electricity_renewables', 'Renewables - % electricity', 'value'};
valIdx = [];
for nm = valCandidates
    j = find(strcmp(have, nm{1}), 1);       % exact match first
    if isempty(j), j = find(strcmpi(have, nm{1}), 1); end
    if ~isempty(j), valIdx = j; break; end
end
if isempty(valIdx)
    % Fallback: pick the rightmost numeric column not in {Entity, Code, Year}
    numMask = varfun(@(x) isnumeric(x), T_owid_all, 'OutputFormat','uniform');
    excl = false(size(have));
    entIdx = find(strcmpi(have,'Entity') | strcmpi(have,'Country'), 1);
    if ~isempty(codeIdx), excl(codeIdx) = true; end
    if ~isempty(yearIdx), excl(yearIdx) = true; end
    if ~isempty(entIdx),  excl(entIdx)  = true; end
    cand = find(numMask & ~excl);
    if ~isempty(cand), valIdx = cand(end); end
end

assert(~isempty(codeIdx) && ~isempty(yearIdx) && ~isempty(valIdx), ...
    "OWID CSV: could not resolve Code/Year/Value columns. Have: %s", strjoin(have, ", "));

% Standardize to {Code, Year, REN_share_pct}
T_owid = T_owid_all(:, [codeIdx, yearIdx, valIdx]);
T_owid.Properties.VariableNames = {'Code','Year','REN_share_pct'};
T_owid.Code = string(T_owid.Code);
T_owid.Year = double(T_owid.Year);
T_owid.REN_share_pct = double(T_owid.REN_share_pct);

% Filter to target countries and clean
T_owid = T_owid(ismember(T_owid.Code, countries) & ...
                ~isnan(T_owid.Year) & ~isnan(T_owid.REN_share_pct), :);
T_owid = sortrows(T_owid, 'Year', 'ascend');

%% ---- Helper handles ----
fetchWB = @(ctyISO3, indicator) wb_fetch_indicator(ctyISO3, indicator, rawDir);

%% ---- 0.2 Fetch WB + merge per country, then build transforms ----
for c = countries
    cISO3 = upper(string(c));
    fprintf('\n---- Country %s ----\n', cISO3);

    % World Bank pulls (each writes WB_<INDICATOR>_<CTY>.csv + META)
    T_gdp = fetchWB(cISO3, WB.GDP);
    T_k   = fetchWB(cISO3, WB.K);
    T_x   = fetchWB(cISO3, WB.X);
    T_m   = fetchWB(cISO3, WB.IMP);

    % OWID subset
    Tr = T_owid(strcmpi(T_owid.Code, cISO3), {'Year','REN_share_pct'});

    % Raw merged (outer join for visibility of gaps)
    T = outerjoin(Tr, T_gdp, 'Keys', 'Year','MergeKeys',true);
    T = outerjoin(T,  T_k,   'Keys', 'Year','MergeKeys',true);
    T = outerjoin(T,  T_x,   'Keys', 'Year','MergeKeys',true);
    T = outerjoin(T,  T_m,   'Keys', 'Year','MergeKeys',true);
    T.Properties.VariableNames = {'Year','REN_share_pct','GDP_2015USD','K_2015USD','X_2015USD','IMP_2015USD'};
    T = sortrows(T, 'Year', 'ascend');

    rawCsv = fullfile(rawDir, sprintf("merged_%s_raw.csv", cISO3));
    writetable(T, rawCsv);
    fprintf("Saved merged raw for %s: %s\n", cISO3, rawCsv);

    % Provenance for merged raw
    fid = fopen(fullfile(rawDir, sprintf("merged_%s_META.txt", cISO3)), 'w');
    fprintf(fid, "Merged raw for %s\n", cISO3);
    fprintf(fid, "Sources:\n- OWID share-electricity-renewables (percent, 0–100)\n");
    fprintf(fid, "- World Bank WDI indicators (constant 2015 US$): %s, %s, %s, %s\n", ...
            WB.GDP, WB.K, WB.X, WB.IMP);
    fprintf(fid, "Window target: %d–%d\nFetched: %s\n", analysisWindow(1), analysisWindow(2), datestr(now,'yyyy-mm-dd HH:MM:SS'));
    fprintf(fid, "Licenses: OWID (CC BY), WDI (CC BY 4.0)\n");
    fclose(fid);

%% ---- 0.5 Build analysis-ready transforms (deterministic) ----
% 1) Keep analysis window and rows with all needed inputs present
T2 = T(T.Year >= analysisWindow(1) & T.Year <= analysisWindow(2), :);
T2 = rmmissing(T2, 'DataVariables', {'REN_share_pct','GDP_2015USD','K_2015USD','X_2015USD','IMP_2015USD'});

% 2) Guards for log transforms: flows must be >0; REN share ∈ (0,1]
assert(all(T2.GDP_2015USD > 0),  'Non-positive GDP values encountered.');
assert(all(T2.K_2015USD   > 0),  'Non-positive K (GFCF) values encountered.');
assert(all(T2.X_2015USD   > 0),  'Non-positive Exports values encountered.');
assert(all(T2.IMP_2015USD > 0),  'Non-positive Imports values encountered.');

eps_floor = 1e-8;                        % guard against log(0) for shares only
s = T2.REN_share_pct / 100;              % percent -> fraction
s = max(s, eps_floor);                   % clamp lower bound only

% 3) Differenced log transforms for VAR in differences
out = table();
out.Year = T2.Year(2:end);
out.gdp  = diff(log(T2.GDP_2015USD));    % Δln GDP
out.k    = diff(log(T2.K_2015USD));      % Δln K (GFCF)
dx       = diff(log(T2.X_2015USD));      % Δln Exports
dm       = diff(log(T2.IMP_2015USD));    % Δln Imports
out.nx   = dx - dm;                      % Δln(X) − Δln(IMP)
out.ren  = diff(log(s));                 % Δln(REN share fraction)

% 4) Write analysis-ready CSV + a transform contract TXT
tCsv = fullfile(rawDir, sprintf("merged_%s_transforms.csv", cISO3));
writetable(out, tCsv);
fprintf("Saved analysis-ready transforms for %s: %s (T=%d)\n", cISO3, tCsv, height(out));

tTxt = fullfile(rawDir, sprintf("merged_%s_TRANSFORMS.txt", cISO3));
fid = fopen(tTxt,'w');
fprintf(fid, "Transforms for %s (analysis window %d–%d)\n", cISO3, analysisWindow(1), analysisWindow(2));
fprintf(fid, "Variables (Δln units unless noted):\n");
fprintf(fid, "  gdp  = Δln(GDP_2015USD)\n");
fprintf(fid, "  k    = Δln(K_2015USD)  [GFCF]\n");
fprintf(fid, "  nx   = Δln(X_2015USD) − Δln(IMP_2015USD)\n");
fprintf(fid, "  ren  = Δln(REN_share), where REN_share = REN_share_pct/100 ∈ (0,1]\n");
nClamp = sum((T2.REN_share_pct/100) < eps_floor);
fprintf(fid, "Notes: rows with missing inputs dropped before differencing; %d REN-share values clamped at %g to avoid log(0).\n", nClamp, eps_floor);
fclose(fid);

end

disp('Step 0 complete.');

%% ----------------------- local functions -----------------------

function T = wb_fetch_indicator(cty_iso3, indicator, rawDir)
% Robust World Bank v2 indicator fetch with pagination.
% Tries ISO-3 then ISO-2; writes CSV + META; returns table(Year, Value).

    cISO3 = upper(string(cty_iso3));
    cISO2 = iso3_to_iso2(cISO3);   % known map for GRC/BEL; falls back to first 2 chars

    candCodes = unique([cISO3, cISO2], 'stable');
    dateWin  = "1960:2025";
    perPage  = 2000;  % bigger than #years but we still paginate robustly
    opts     = weboptions('Timeout', 60, 'ContentType', 'json');

    years = []; vals = []; got = false; lastUrl = "";

    for cc = candCodes
        page = 1; localYears = []; localVals = [];
        while true
            url = sprintf("https://api.worldbank.org/v2/country/%s/indicator/%s?format=json&date=%s&per_page=%d&page=%d", ...
                          lower(cc), upper(indicator), dateWin, perPage, page);
            lastUrl = url;
            raw = webread(url, opts);
            [meta, data] = parse_wb_payload(raw);
            [y,v] = extract_wb_year_values(data);
            localYears = [localYears; y]; localVals = [localVals; v];

            pages = 1;
            if isstruct(meta) && isfield(meta,'pages') && ~isempty(meta.pages)
                pages = max(1, double(meta.pages));
            end
            if page >= pages, break; end
            page = page + 1;
        end
        if ~isempty(localYears)
            years = localYears; vals = localVals; got = true; break;
        end
    end

    if ~got
        error("WB fetch failed for %s/%s. Last URL: %s", cISO3, indicator, lastUrl);
    end

    ok = ~isnan(years) & ~isnan(vals);
    T  = table(years(ok), vals(ok), 'VariableNames', {'Year','Value'});
    T  = sortrows(unique(T, 'rows'), 'Year');

    outCsv = fullfile(rawDir, sprintf("WB_%s_%s.csv", upper(indicator), cISO3));
    writetable(T, outCsv);

    % Provenance TXT
    metaTxt = strrep(outCsv, '.csv', '_META.txt');
    fid = fopen(metaTxt,'w');
    fprintf(fid, "World Bank WDI %s | country candidates=%s | window=%s | per_page=%d | fetched=%s\n", ...
            upper(indicator), strjoin(cellstr(candCodes),','), dateWin, perPage, datestr(now,'yyyy-mm-dd HH:MM:SS'));
    fprintf(fid, "License: CC BY 4.0\n");
    fclose(fid);
end

function [meta, data] = parse_wb_payload(raw)
    if iscell(raw) && numel(raw) >= 2
        meta = raw{1}; data = raw{2};
        if isstruct(meta) && isfield(meta,'message')
            error("WB API error: %s", jsonencode(meta));
        end
    else
        error("Unexpected WB payload type: %s", class(raw));
    end
end

function [years, vals] = extract_wb_year_values(data)
    years = []; vals = [];
    if isstruct(data)
        for k = 1:numel(data)
            if ~isfield(data(k),'date') || ~isfield(data(k),'value'), continue; end
            y = str2double(data(k).date); v = data(k).value;
            if ischar(v) || isstring(v), v = str2double(v); end
            if ~isnan(y) && ~isempty(v) && ~isnan(v)
                years(end+1,1) = y; vals(end+1,1)  = double(v);
            end
        end
    elseif iscell(data)
        for k = 1:numel(data)
            rec = data{k};
            if ~isstruct(rec) || ~isfield(rec,'date') || ~isfield(rec,'value'), continue; end
            y = str2double(rec.date); v = rec.value;
            if ischar(v) || isstring(v), v = str2double(v); end
            if ~isnan(y) && ~isempty(v) && ~isnan(v)
                years(end+1,1) = y; vals(end+1,1)  = double(v);
            end
        end
    else
        error("WB data block type: %s", class(data));
    end
end

function iso2 = iso3_to_iso2(iso3)
% Minimal ISO-3 → ISO-2 map for our countries; fallback = first two letters.
    iso3 = upper(string(iso3));
    switch iso3
        case "GRC", iso2 = "GR";
        case "BEL", iso2 = "BE";
        otherwise,  iso2 = extractBefore(iso3, 3);  % crude fallback
    end
end