Fix errors and bugs in the given code, without changing its intended functionality. ONLY return the fixed code and nothing else. Enclose the entire code in a single code block.



function [best_theta_est, best_logL, elapsed_time, H_t, p_values, coverage] = VTNLscalarBEKK_V2(X, initial_params, non_linear)
    % Start timer
    tic;

    % Define parameter grid based on model type (with or without non-linear component)
    if non_linear
        alpha_grid = linspace(0.001, 0.2, 3);
        beta_grid = linspace(0.5, 0.99, 3);
        phi_grid = linspace(-0.1, 0.1, 3);
        gamma_grid = linspace(0.001, 50, 3);
        c_grid = linspace(-0.09, 0.09, 3);

        % Create parameter combinations
        param_combinations = combvec(alpha_grid, beta_grid, phi_grid, gamma_grid, c_grid)';
    else
        alpha_grid = linspace(0.001, 0.2, 3);
        beta_grid = linspace(0.5, 0.99, 3);

        % Create parameter combinations
        param_combinations = combvec(alpha_grid, beta_grid)';
    end

    num_combinations = size(param_combinations, 1);

    % If initial_params is provided, add it to the parameter combinations
    if nargin > 1 && ~isempty(initial_params)
        param_combinations = [param_combinations; initial_params];
        num_combinations = num_combinations + 1;
    end

    % Initialize results storage
    logL_results = -Inf(num_combinations, 1);
    theta_results = zeros(num_combinations, size(param_combinations, 2));

    % Use parallel for loop for computation
    parfor i = 1:num_combinations
        theta_init = param_combinations(i, :);

        % Transform initial parameters to ensure validity during optimization
        theta_init_transformed = transform_parameters(theta_init);

        % Perform optimization with fminunc
        options = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'off');
        try
            [theta_est_transformed, fval] = fminunc(@(theta) -log_likelihood_vtnlsbekk(X, inverse_transform_parameters(theta), non_linear), theta_init_transformed, options);
            theta_est = inverse_transform_parameters(theta_est_transformed);
        catch
            % Skip if optimization fails
            continue;
        end

        % Store results
        logL_results(i) = -fval;
        theta_results(i, :) = theta_est;
    end

    % Find best result
    [best_logL, best_idx] = max(logL_results);
    best_theta_est = theta_results(best_idx, :);

    % End timer and calculate elapsed time
    elapsed_time = toc;
    elapsed_time = elapsed_time / 60;

    % Calculate H_t for the best parameters
    [~, H_t] = log_likelihood_vtnlsbekk(X, best_theta_est, non_linear);

    % Calculate p-values and coverage
    [p_values, coverage] = calculate_p_values_and_coverage(X, best_theta_est, non_linear);
end

function theta_transformed = transform_parameters(theta)
    % Example transformation: log transformation to ensure positivity
    theta_transformed = log(theta);
end

function theta = inverse_transform_parameters(theta_transformed)
    % Inverse of the log transformation
    theta = exp(theta_transformed);
end

function [logL, H_t] = log_likelihood_vtnlsbekk(X, theta, non_linear)
    % Unpack parameters
    alpha = theta(1);
    beta = theta(2);

    if non_linear
        phi = theta(3);
        gamma = theta(4);
        c = theta(5);
    else
        phi = 0;
        gamma = 0;
        c = 0;
    end

    % Ensure parameters meet conditions
    if beta <= 0 || alpha <= 0 || alpha + phi <= 0 || alpha + phi + beta >= 1 || gamma < 0
        logL = -Inf;
        return;
    end

    % Dimensions
    [T, ~] = size(X);
    H = cov(X);  % Initial unconditional covariance matrix
    H_t = repmat(H, [1, 1, T]);  % Initial conditional covariance matrix

    logL = 0;

    % Compute log-likelihood
    for t = 2:T
        if non_linear
            G_tc = 1 / (1 + exp(-gamma * (t / T - c))); % Calculate G_tc normally
            H_t(:, :, t) = (1 - alpha - phi * ((T + 1) / (2 * T)) - beta) * H ...
                            + (alpha + phi * G_tc) * (X(t-1, :)' * X(t-1, :)) ...
                            + beta * H_t(:, :, t-1);
        else
            H_t(:, :, t) = (1 - alpha - beta) * H ...
                            + alpha * (X(t-1, :)' * X(t-1, :)) ...
                            + beta * H_t(:, :, t-1);
        end

        % Check if covariance matrix is invertible
        if ~isfinite(det(H_t(:, :, t)))
            logL = -Inf;
            return;
        end

        logL = logL + log(det(H_t(:, :, t))) + X(t, :) * (H_t(:, :, t)  X(t, :)');
    end

    logL = -0.5 * logL;
end

function [p_values, coverage] = calculate_p_values_and_coverage(X, theta_est, non_linear)
    % Calculate standard errors, p-values, and 95% coverage for alpha and beta
    [~, H_t] = log_likelihood_vtnlsbekk(X, theta_est, non_linear);
    T = size(X, 1);

    % Example calculation, you might need to adjust this based on your actual implementation
    alpha = theta_est(1);
    beta = theta_est(2);

    % Assuming normal distribution for p-value calculation
    alpha_se = std(squeeze(H_t(1,1,:))) / sqrt(T);
    beta_se = std(squeeze(H_t(2,2,:))) / sqrt(T);

    % Calculate z-scores
    alpha_z = alpha / alpha_se;
    beta_z = beta / beta_se;

    % Calculate p-values
    alpha_p = 2 * (1 - normcdf(abs(alpha_z)));
    beta_p = 2 * (1 - normcdf(abs(beta_z)));

    p_values = [alpha_p, beta_p];

    % Calculate 95% coverage
    alpha_coverage = (alpha - 1.96 * alpha_se < alpha) && (alpha + 1.96 * alpha_se > alpha);
    beta_coverage = (beta - 1.96 * beta_se < beta) && (beta + 1.96 * beta_se > beta);

    coverage = [alpha_coverage, beta_coverage];