function [data_moments, m_data] = get_data_moments(data, obs_var, inv_order_var, matched_moments_)
    mom_nbr = size(matched_moments_, 1);
    T = size(data, 1);
    
    % Extract variables (match varobs order)
    c_data = data(:, 1);
    invest_data = data(:, 2);
    log_y_data = data(:, 3);
    nx_data = data(:, 4);
    delta_y_data = data(:, 5);
    
    % Demean variables
    c_dm = c_data - mean(c_data);
    i_dm = invest_data - mean(invest_data);
    y_dm = log_y_data - mean(log_y_data);
    nx_dm = nx_data - mean(nx_data);
    dy_dm = delta_y_data - mean(delta_y_data);
    
    % Compute standard deviations
    std_y = std(log_y_data);
    std_c = std(c_data);
    std_i = std(invest_data);
    std_nx = std(nx_data);
    std_dy = std(delta_y_data);
    var_y = std_y^2;
    var_c = std_c^2;
    var_i = std_i^2;
    var_nx = std_nx^2;
    var_dy = std_dy^2;
    %{
    % Compute ALL 11 moments - ORIGINAL
    all_data_moments = zeros(11, 1);
    all_data_moments(1) = std_y * 100;                  % σ(y) in %
    all_data_moments(2) = std_dy * 100;                 % σ(Δy) in %
    all_data_moments(3) = std_c / std_y;                % σ(c)/σ(y)
    all_data_moments(4) = std_i / std_y;                % σ(I)/σ(y)
    all_data_moments(5) = std_nx / std_y;               % σ(nx)/σ(y)
    all_data_moments(6) = corr(log_y_data(2:end), log_y_data(1:end-1));
    all_data_moments(7) = corr(delta_y_data(2:end), delta_y_data(1:end-1));
    all_data_moments(8) = corr(log_y_data, nx_data);
    all_data_moments(9) = corr(log_y_data, c_data);
    all_data_moments(10) = corr(log_y_data, invest_data);
    all_data_moments(11) = mean(delta_y_data);
    
    % Create m_data for HAC-robust variance estimation
    all_m_data = zeros(T, 11);
    
     % Moment 1: σ(y)
    all_m_data(:,1) = ((y_dm.^2 - var_y) ./ (2 * std_y)) * 100;
    
    % Moment 2: σ(Δy)
    all_m_data(:,2) = ((dy_dm.^2 - var_dy) ./ (2 * std_dy)) * 100;
    
    % Moment 3: σ(c)/σ(y)
    all_m_data(:,3) = ...
        ( (c_dm.^2 - var_c) ./ (2 * std_c * std_y) ) ...
      - ( (std_c ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    % Moment 4: σ(I)/σ(y)
    all_m_data(:,4) = ...
        ( (i_dm.^2 - var_i) ./ (2 * std_i * std_y) ) ...
      - ( (std_i ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    % Moment 5: σ(nx)/σ(y)
    all_m_data(:,5) = ...
        ( (nx_dm.^2 - var_nx) ./ (2 * std_nx * std_y) ) ...
      - ( (std_nx ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    
    % Moment 6: ρ(y,y₋₁) - correlation
    all_m_data(2:end, 6) = (y_dm(1:end-1) .* y_dm(2:end)) / (std_y^2) - all_data_moments(6);
    all_m_data(1, 6) = mean(all_m_data(2:end, 6));
    
    % Moment 7: ρ(Δy,Δy₋₁)
    all_m_data(2:end, 7) = (dy_dm(1:end-1) .* dy_dm(2:end)) / (std_dy^2) - all_data_moments(7);
    all_m_data(1, 7) = mean(all_m_data(2:end, 7));
    
    % Moment 8: ρ(y,nx)
    all_m_data(:, 8) = (y_dm .* nx_dm) / (std_y * std_nx) - all_data_moments(8);
    
    % Moment 9: ρ(y,c)
    all_m_data(:, 9) = (y_dm .* c_dm) / (std_y * std_c) - all_data_moments(9);
    
    % Moment 10: ρ(y,I)
    all_m_data(:, 10) = (y_dm .* i_dm) / (std_y * std_i) - all_data_moments(10);
    
    % Moment 11: E[Δy] - simple mean
    all_m_data(:, 11) = delta_y_data - all_data_moments(11);

    % ============================================================
    % RETURN THE CORRECT MOMENT SET
    % ============================================================
    if mom_nbr >= 11
        data_moments = all_data_moments;
        m_data = all_m_data;
    elseif mom_nbr == 1
        data_moments = all_data_moments(11);
        m_data = all_m_data(:, 11);
    else
        data_moments = all_data_moments(1:10);
        m_data = all_m_data(:, 1:10);
    end
    %}
    
    
    all_data_moments = zeros(10, 1);
    all_data_moments(1) = std_y * 100;                  % σ(y)
    all_data_moments(2) = std_dy * 100;                 % σ(Δy)
    all_data_moments(3) = std_c / std_y;                % σ(c)/σ(y)
    all_data_moments(4) = std_i / std_y;                % σ(I)/σ(y)
    all_data_moments(5) = std_nx / std_y;               % σ(nx)/σ(y)
    all_data_moments(6) = corr(log_y_data(2:end), log_y_data(1:end-1));
    all_data_moments(7) = corr(delta_y_data(2:end), delta_y_data(1:end-1));
    all_data_moments(8) = corr(log_y_data, nx_data);
    all_data_moments(9) = corr(log_y_data, c_data);
    all_data_moments(10) = corr(log_y_data, invest_data);
    
    % ============================================================
    % COMPUTE m_data FOR HAC VARIANCE
    % ============================================================
    all_m_data = zeros(T, 10);
    
    % Moment 1: σ(y)
    all_m_data(:,1) = ((y_dm.^2 - var_y) ./ (2 * std_y)) * 100;
    
    % Moment 2: σ(Δy)
    all_m_data(:,2) = ((dy_dm.^2 - var_dy) ./ (2 * std_dy)) * 100;
    
    % Moment 3: σ(c)/σ(y)
    all_m_data(:,3) = ...
        ( (c_dm.^2 - var_c) ./ (2 * std_c * std_y) ) ...
      - ( (std_c ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    % Moment 4: σ(I)/σ(y)
    all_m_data(:,4) = ...
        ( (i_dm.^2 - var_i) ./ (2 * std_i * std_y) ) ...
      - ( (std_i ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    % Moment 5: σ(nx)/σ(y)
    all_m_data(:,5) = ...
        ( (nx_dm.^2 - var_nx) ./ (2 * std_nx * std_y) ) ...
      - ( (std_nx ./ (2 * std_y^3)) .* (y_dm.^2 - var_y) );
    
    % Moment 6: ρ(y,y₋₁)
    all_m_data(2:end,6) = (y_dm(1:end-1) .* y_dm(2:end)) / (std_y^2) - all_data_moments(6);
    all_m_data(1,6) = mean(all_m_data(2:end,6));
    
    % Moment 7: ρ(Δy,Δy₋₁)
    all_m_data(2:end,7) = (dy_dm(1:end-1) .* dy_dm(2:end)) / (std_dy^2) - all_data_moments(7);
    all_m_data(1,7) = mean(all_m_data(2:end,7));
    
    % Moment 8: ρ(y,nx)
    all_m_data(:,8) = (y_dm .* nx_dm) / (std_y * std_nx) - all_data_moments(8);
    
    % Moment 9: ρ(y,c)
    all_m_data(:,9) = (y_dm .* c_dm) / (std_y * std_c) - all_data_moments(9);
    
    % Moment 10: ρ(y,I)
    all_m_data(:,10) = (y_dm .* i_dm) / (std_y * std_i) - all_data_moments(10);
    
    % 🚫 DELETE OR COMMENT OUT THIS OLD PART
    %{
    % Moment 11: E[Δy]
    % all_m_data(:,11) = delta_y_data - all_data_moments(11);
    %}
    
    % ============================================================
    % RETURN CORRECT MOMENT SET (SIMPLE AND SAFE)
    % ============================================================
    data_moments = all_data_moments(1:mom_nbr);
    m_data       = all_m_data(:, 1:mom_nbr);
    %{
    % ============================================
    % Compute ALL 4 data moments
    % ============================================
    all_data_moments = zeros(3, 1);
    all_data_moments(1) = std_y * 100;                  % σ(y) in %
    all_data_moments(2) = std_c / std_y;                % σ(c)/σ(y)
    
    all_data_moments(3) = corr(log_y_data, c_data);
    m_data = eye(3);
    
    % Determine which moments are requested based on matched_moments_
    % Check if this is the "first moments" call or "other moments" call
    if mom_nbr == 3
        % Full call - return all
        data_moments = all_data_moments;
    elseif mom_nbr == 1
        % Only first moment (mean(delta_y))
        data_moments = all_data_moments(3);
    else
        % Other moments (the remaining 10)
        data_moments = all_data_moments(1:2);
    end
    %}


end