function KODE_FINAL
    clc; clear; close all;
    
    % --- 1. DATA IMPORT & PROCESSING ---
    FileName = 'XYZ_Stock_Price_Data.xlsx'; 
    try
        raw_data = readmatrix(FileName);
        Prices = raw_data(:, 5);
        Prices = Prices(~isnan(Prices));
        Prices = flipud(Prices);
        
        % Return & Volatility Calculation
        Returns = diff(log(Prices));
        Mean_Ret = mean(Returns);
        Std_Ret  = std(Returns);
        Vol_Ann  = Std_Ret * sqrt(252);
        
        % --- 2. PARAMETER DEFINITIONS ---
        S0          = 189.11;       
        r           = 0.0352;       
        sigma       = 0.4991;       
        D_amount    = 0.01;         
        Days_Expiry = 72;           
        Days_Div    = 64;           
        
        T_Exp = Days_Expiry/365;
        T_Div = Days_Div/365;
        
        % Dividend Correction
        S_star = S0 - D_amount * exp(-r * T_Div);
        
        % Display Parameters
        fprintf('Underlying Asset    : XYZ Tech Corp (XYZ)\n');
        fprintf('Spot Price (S0)     : %.2f USD\n', S0);
        fprintf('Adjusted Price (S*) : %.5f USD\n', S_star);
        fprintf('Risk-Free Rate (r)  : %.2f%%\n', r*100);
        fprintf('Volatility (sigma)  : %.4f\n\n', sigma);
        
        % --- 3. GRAPHICAL OUTPUT ---
        
        % Figure 1: Historical Prices
        figure('Color','w');
        plot(Prices, 'b-', 'LineWidth', 1.5);
        title('Historical Stock Price Movement (XYZ Tech Corp)');
        subtitle('Observation Period: 1 Year (Daily Closing Prices)');
        xlabel('Trading Days'); ylabel('Stock Price (USD)');
        grid on; xlim([1 length(Prices)]);
        
        % Figure 2: Normality Test 
        figure('Color','w');
        histogram(Returns, 40, 'Normalization', 'pdf', 'FaceColor', [0.7 0.7 0.7], 'EdgeColor', 'w'); 
        hold on;
        x_val = linspace(min(Returns), max(Returns), 100);
        y_val = (1 ./ (Std_Ret * sqrt(2*pi))) .* exp(-0.5 * ((x_val - Mean_Ret) ./ Std_Ret).^2);
        plot(x_val, y_val, 'r-', 'LineWidth', 2.5);
        title('Log-Return Distribution vs Theoretical Normal Curve');
        xlabel('Daily Log Return'); ylabel('Probability Density');
        legend({'Empirical Data', 'Theoretical Normal Dist.'}, 'Location', 'NorthWest');
        grid on;
        
        % --- 4. CONVERGENCE TEST ---
        K_test = 190; 
        [BS_Bench, ~] = BlackScholes(S_star, K_test, T_Exp, r, sigma);
        N_values = [10, 20, 50, 100, 200, 500, 1000];
        
        fprintf('N        Numerical Price (USD)   Abs. Error\n');
        fprintf('-------------------------------------------\n');
        for i = 1:length(N_values)
            n = N_values(i);
            P_Simp = Simpson13(S_star, K_test, T_Exp, r, sigma, 'call', n, 4);
            err = abs(P_Simp - BS_Bench);
            fprintf('%-8d %-23.5f %-15.8f\n', n, P_Simp, err);
        end
        fprintf('\n');

        % --- 5. VALUATION TABLES ---
        Strikes  = [140; 155; 190; 220; 240];
        Last_Price_Call = [51.40; 38.85; 14.67; 4.60; 1.89];
        Last_Price_Put  = [1.40; 3.11; 14.00; 34.10; 51.40];
        
        Results_Call = {}; 
        Results_Put = {};
        Total_APE_Call = 0; 
        Total_APE_Put = 0;
        
        for i = 1:length(Strikes)
            K = Strikes(i);
            
            % Determine Moneyness
            Moneyness_Call = GetMoneyness(S0, K, 'call');
            Moneyness_Put  = GetMoneyness(S0, K, 'put');
            
            % Calculations
            [bs_c, bs_p] = BlackScholes(S_star, K, T_Exp, r, sigma);
            simp_c = Simpson13(S_star, K, T_Exp, r, sigma, 'call', 200, 4);
            simp_p = Simpson13(S_star, K, T_Exp, r, sigma, 'put', 200, 4);
            
            % Error Calculation (APE)
            ape_c  = abs(Last_Price_Call(i) - simp_c) / Last_Price_Call(i) * 100;
            ape_p  = abs(Last_Price_Put(i) - simp_p) / Last_Price_Put(i) * 100;
            
            Total_APE_Call = Total_APE_Call + ape_c;
            Total_APE_Put  = Total_APE_Put + ape_p;
            
            % Store Data
            Results_Call(i,:) = {Moneyness_Call, K, Last_Price_Call(i), simp_c, bs_c, ape_c};
            Results_Put(i,:)  = {Moneyness_Put,  K, Last_Price_Put(i),  simp_p, bs_p, ape_p};
        end
        
        MAPE_Call = Total_APE_Call / length(Strikes);
        MAPE_Put  = Total_APE_Put  / length(Strikes);
        
        PrintTable('CALL OPTION VALUATION', Results_Call, MAPE_Call);
        PrintTable('PUT OPTION VALUATION',  Results_Put,  MAPE_Put);
        
    catch ME
        fprintf('Error: %s\n', ME.message);
    end
end

% --- HELPER FUNCTIONS ---

function PrintTable(Title, Data, MAPE)
    fprintf('%s\n', Title);
    fprintf('%-10s %-18s %-16s %-28s %-28s %-8s\n', ...
            'Zone', 'Strike Price (USD)', 'Last Price (USD)', ...
            'Numerical Option Price (USD)', 'Analytical Option Price (USD)', 'APE (%)');
    
    fprintf('------------------------------------------------------------------------------------------------------------------------\n');
    
    for i = 1:size(Data,1)
        fprintf('%-10s %-18d %-16.2f %-28.5f %-28.5f %-8.2f\n', ...
                Data{i,1}, Data{i,2}, Data{i,3}, Data{i,4}, Data{i,5}, Data{i,6});
    end
    fprintf('------------------------------------------------------------------------------------------------------------------------\n');
    fprintf('Mean APE: %.2f%%\n\n', MAPE);
end

function zone = GetMoneyness(S, K, type)
    ratio = S/K;
    if ratio > 1.01 
        if strcmpi(type,'call'), zone='ITM'; else, zone='OTM'; end
    elseif ratio < 0.99 
        if strcmpi(type,'call'), zone='OTM'; else, zone='ITM'; end
    else 
        zone = 'ATM';
    end
end

function [call, put] = BlackScholes(S, K, T, r, sigma)
    d1 = (log(S/K) + (r + 0.5*sigma^2)*T) / (sigma*sqrt(T));
    d2 = d1 - sigma*sqrt(T);
    call = S * MyNorm(d1) - K * exp(-r*T) * MyNorm(d2);
    put  = K * exp(-r*T) * MyNorm(-d2) - S * MyNorm(-d1);
end

function val = Simpson13(S, K, T, r, sigma, type, N, limit)
    mu = (r - 0.5*sigma^2)*T;
    if strcmpi(type,'call')
        a = K; b = S*limit;
    else
        a = 1e-5; b = K;
    end
    h = (b-a)/N; 
    sum_odd = 0; sum_even = 0;
    for i=1:N-1
        x = a + i*h; 
        y = Integrand(x,S,K,mu,sigma,T,type);
        if mod(i,2)==1, sum_odd = sum_odd + y; else, sum_even = sum_even + y; end
    end
    integral_val = (h/3) * (Integrand(a,S,K,mu,sigma,T,type) + 4*sum_odd + 2*sum_even + Integrand(b,S,K,mu,sigma,T,type));
    val = exp(-r*T) * integral_val;
end

function out = Integrand(St, S, K, mu, sigma, T, type)
    if St <= 0, out = 0; return; end
    pdf = (1/(St*sigma*sqrt(2*pi*T))) * exp(-((log(St/S)-mu)^2)/(2*sigma^2*T));
    if strcmpi(type,'call'), pay = max(St - K, 0); else, pay = max(K - St, 0); end
    out = pay * pdf;
end

function p = MyNorm(x)
    p = 0.5 * (1 + erf(x/sqrt(2))); 
end