% 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
