% Joint PDF of Area Load (L) and Transfer Limit (TL)

% --------------------------------------------------

% INPUTS -------------------------------------------------

% L – column vector (N×1) of historical area-load values

% TL – column vector (N×1) of the corresponding transfer-limit values

%

% OUTPUTS ------------------------------------------------

% X1, X2 – evaluation grid (load, TL) for plotting / lookup

% fJoint – matrix of joint-pdf estimates at each (X1(i,j), X2(i,j))

%% 1. Load or assign your data ------------------------------------------

% Replace these with your actual vectors

L = load('areaLoad_MW.mat' ,'-mat'); % e.g. struct with field L

TL = load('transferLimit_MW.mat','-mat'); % struct with field TL

L = L.L(:); % ensure column shape

TL = TL.TL(:);

data = [L TL]; % N×2 matrix for ksdensity

%% 2. Build an evaluation grid ------------------------------------------

nGrid = 150; % resolution of the grid

x1 = linspace(min(L), max(L), nGrid); % load-axis points

x2 = linspace(min(TL), max(TL), nGrid); % TL-axis points

[X1,X2] = meshgrid(x1,x2);

gridPts = [X1(:) X2(:)]; % flatten for ksdensity

%% 3. 2-D kernel density estimate of the joint PDF -----------------------

% ‘ksdensity’ uses a product Gaussian kernel; bandwidth is

% automatically selected (Silverman’s rule) unless you override it.

f = ksdensity(data, gridPts, ...

          'Function',  'pdf', ...

          'Support',   'unbounded');    % returns vector length nGrid^2

fJoint = reshape(f, nGrid, nGrid); % back to 2-D matrix

%% 4. (Optional) Plot the surface ----------------------------------------

figure;

surf(X1, X2, fJoint, 'EdgeColor', 'none');

xlabel('Area Load (MW)');

ylabel('Transfer Limit (MW)');

zlabel('Joint PDF f_{L,TL}(l, tl)');

title('Kernel Joint-PDF of Transfer Limit vs. Load');

view([30 40]); % nice 3-D viewing angle

colormap parula