function R = weightedcorrs(Y, w)
%
% WEIGHTEDCORRS returns a symmetric matrix R of weighted correlation
% coefficients calculated from an input T-by-N matrix Y whose rows are
% observations and whose columns are variables and an input T-by-1 vector
% w of weights for the observations. This function may be a valid
% alternative to CORRCOEF if observations are not all equally relevant
% and need to be weighted according to some theoretical hypothesis or
% knowledge.
%
% R = WEIGHTEDCORRS(Y, w) returns a positive semidefinite matrix R,
% i.e. all its eigenvalues are non-negative (see Example 1).
%
% WEIGHTEDCORRS is such that
% WEIGHTEDCORRS(Y, w) == WEIGHTEDCORRS(a * Y + b, w)
% where a and b are two real numbers (see Example 1).
%
% Furthermore, the result provided by the function doesn't change if the
% unit system of each column of Y is changed through an arbitrary affine
% transformation y = a * x + b, where a and b are two real numbers, with
% a > 0 (see Example 2).
%
% If w = ones(size(Y, 1), 1), no difference exists between
% WEIGHTEDCORRS(Y, w) and CORRCOEF(Y) (see Example 4).
%
%
% % ======================================================================
% % EXAMPLE 0: some common shapes in 2D.
% % ======================================================================
%
% T = 1000; % number of observations
% % CHOOSE WEIGHTS
% alpha = 2 / T;
% w0 = 1 / sum(exp(((1:T) - T) * alpha));
% w = w0 * exp(((1:T) - T) * alpha); % weights: exponential decay
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % a) A line is always a line!
% Y(:, 1) = 1:T;
% Y(:, 2) = rand * Y(:, 1) + rand; % Linear relation
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % b) An horizontal line has a correlation equal to 0
% Y(1:T, 1) = rand;
% Y(:, 2) = 1:T; % Linear relation
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % c) A vertical line has a correlation equal to 0
% Y(:, 1) = 1:T;
% Y(:, 2) = rand; % Linear relation
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % d) A symmetric parabolic shape ... makes a huge difference!
% Y(:, 1) = (1:T) - T / 2 - 0.5;
% Y(:, 2) = rand * Y(:, 1) .^ 2 + rand; % Parabolic relation (symmetric)
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % e) A circle
% angle = (2 * rand(T, 1) - 1) * pi;
% rho = rand(T, 1);
% Y(:, 1) = rho .* cos(angle);
% Y(:, 2) = rho .* sin(angle); % Circle
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % f) An exponential curve
% Y(:, 1) = 1:T;
% Y(:, 2) = exp(3 * (1:T) / T); % Exponential relation
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % g) A logarithm
% Y(:, 1) = 1:T;
% Y(:, 2) = log(1:T); % Logarithmic relation
% r1 = corrcoef(Y);
% r1 = r1(2) % Traditional Correlation
% r2 = weightedcorrs(Y, w);
% r2 = r2(2) % Weighted Correlation
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% % ======================================================================
% % EXAMPLE 1: verify some of the properties for weighted correlations.
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 500; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
%
% % CHOOSE WEIGHTS
% alpha = 2 / T;
% w0 = 1 / sum(exp(((1:T) - T) * alpha));
% w = w0 * exp(((1:T) - T) * alpha); % weights: exponential decay
%
% % COMPUTE CORRELATION MATRIX
% r1 = weightedcorrs(Y, w); % Weighted Correlation matrix for Y
%
% % COMPUTE CORRELATION MATRIX FOR MODIFIED DATA
% a = 1000 * (2 * rand - 1); % Multiplicative coefficient
% b = 1000 * (2 * rand - 1); % Shift
% r2 = weightedcorrs(a * Y + b, w); % Weighted Correlation matrix for modified Y
%
% % FIND RELEVANT INDEXES AND PLOT A SCATTER DIAGRAM
% I = ones(N);
% indexes = find(tril(I, -1)); % Indexes of lower triangular square matrix
%
% figure('Position', [50 50 950 630]);
% plot([-1:1], [-1:1], 'r', 'LineWidth', 10); % 45 Line
% hold on;
% plot(r1(indexes), r2(indexes), '*', 'MarkerSize', 6); % Comparison between the coefficients of the two Matrices
% title('Identical correlations for Y and a * Y + b', ... % title label for the plot
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Correlations for Y', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Correlations for modified Y = a * Y + b', ...
% 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % PLOT HISTOGRAM OF DIFFERENCES BETWEEN COEFFICIENTS
% % OF THE TWO CORRELATION MATRICES
% figure('Position', [50 50 950 630]);
% hist(r1(indexes) - r2(indexes), 100); % Histogram of the differences (which are always negligible)
% title('Identical correlations for Y and a * Y + b', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Negligible differences', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Absolute Frequencies', 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % VERIFY SYMMETRY
% all(all(r1 == r1'))
% all(all(r2 == r2'))
%
% % VERIFY POSITIVENESS FOR ALL EIGENVALUES
% all(eig(r1) >= 0)
% all(eig(r2) >= 0)
%
%
% % ======================================================================
% % EXAMPLE 2: the unit system of each column of Y is changed through
% % arbitrary affine transformations.
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 500; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
%
% % CHOOSE WEIGHTS
% alpha = 2 / T;
% w0 = 1 / sum(exp(((1:T) - T) * alpha));
% w = w0 * exp(((1:T) - T) * alpha); % weights: exponential decay
%
% % COMPUTE CORRELATION MATRIX
% r1 = weightedcorrs(Y, w); % Weighted Correlation matrix for Y
%
% % COMPUTE CORRELATION MATRIX FOR MODIFIED DATA
% a = 1000 * rand(1, N); % Multiplicative coefficients (positive)
% b = 1000 * (2 * rand(1, N) - 1); % Shifts
% r2 = weightedcorrs(repmat(a, T, 1) .* Y + repmat(b, T, 1), w); % Weighted Correlation matrix for modified Y
%
% % FIND RELEVANT INDEXES AND PLOT A SCATTER DIAGRAM
% I = ones(N);
% indexes = find(tril(I, -1)); % Indexes of lower triangular square matrix
%
% figure('Position', [50 50 950 630]);
% plot([-1:1], [-1:1], 'r', 'LineWidth', 10); % 45 Line
% hold on;
% plot(r1(indexes), r2(indexes), '*', 'MarkerSize', 6); % Comparison between the coefficients of the two Matrices
% str = 'Identical correlations after arbitrary affine transformations';
% title(str, 'FontSize', 20, 'FontWeight', 'Bold'); % title label for the plot
% xlabel('Correlations for Y', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Correlations for modified Y', ...
% 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % PLOT HISTOGRAM OF DIFFERENCES BETWEEN COEFFICIENTS
% % OF THE TWO CORRELATION MATRICES
% figure('Position', [50 50 950 630]);
% hist(r1(indexes) - r2(indexes), 100); % Histogram of the differences (which are always negligible)
% title('Identical correlations for Y and modified Y', ... % title label for the plot
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Negligible differences', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Absolute Frequencies', 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
%
% % ======================================================================
% % Example 3: differences with respect to traditional correlations.
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 500; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
%
% % CHOOSE WEIGHTS AND PLOT THEM
% alpha = 2 / T;
% w0 = 1 / sum(exp(((1:T) - T) * alpha));
% w = w0 * exp(((1:T) - T) * alpha); % weights: exponential decay
%
% figure('Position', [50 50 950 630]);
% plot([floor(-T / 2):ceil(T + T / 2)], ... % x values ...
% w0 * exp(((floor(-T / 2):ceil(T + T / 2)) - T) * alpha), ... % ... y values ...
% '.r', 'MarkerSize', 7); % ... plot properties
% hold on;
% plot(w, '.', 'MarkerSize', 30);
% title('Exponential weights assigned to observations', ... % title label for the plot
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Observations', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Weights', 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % COMPUTE CORRELATION MATRICES
% r1 = weightedcorrs(Y, w); % Weighted Correlation matrix
% r2 = corrcoef(Y); % Traditional Correlation matrix
%
% % FIND RELEVANT INDEXES AND PLOT A SCATTER DIAGRAM
% I = ones(N);
% indexes = find(tril(I, -1)); % Indexes of lower triangular square matrix
%
% figure('Position', [50 50 950 630]);
% plot(r1(indexes), r2(indexes), '.'); % Comparison with the Traditional Correlation matrix
% hold on;
% plot([-1:1], [-1:1], 'r', 'LineWidth', 5); % 45 Line
% title('Scatter Diagram for Traditional and Weighted Correlations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Weighted Correlation coefficients', ... % x label ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% ylabel('Traditional Correlation coefficients', ... % y label ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % PLOT HISTOGRAM OF DIFFERENCES BETWEEN COEFFICIENTS
% % OF THE TWO CORRELATION MATRICES
% figure('Position', [50 50 950 630]);
% hist(r2(indexes) - r1(indexes), 100); % Differences between the two correlation matrices
% title('Differences between Traditional and Weighted Correlations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% ylabel('Absolute Frequencies', 'FontSize', 16, 'FontWeight', 'Bold');
% xlabel('Differences between Traditional and Weighted Correlations', ... % x label ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
%
% % ======================================================================
% % Example 4: If weights are uniform, then the traditional correlation
% % matrix is obtained.
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 500; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
%
% % CHOOSE WEIGHTS AND PLOT THEM
% w = ones(T, 1) / T; % weights: uniform
%
% figure('Position', [50 50 950 630]);
% plot(w, '*', 'MarkerSize', 7);
% title('Uniform Weights assigned to observations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Observations', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Weights', 'FontSize', 16, 'FontWeight', 'Bold');
% axis([0 T 0 2 * w(1)]);
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % COMPUTE CORRELATION MATRICES
% r1 = weightedcorrs(Y, w); % Weighted Correlation matrix
% r2 = corrcoef(Y); % Traditional Correlation matrix
%
% % FIND RELEVANT INDEXES AND PLOT A SCATTER DIAGRAM
% I = ones(N);
% indexes = find(tril(I, -1)); % Indexes of lower triangular square matrix
%
% figure('Position', [50 50 950 630]);
% plot([-1:1], [-1:1], 'r', 'LineWidth', 10); % 45 Line
% hold on;
% plot(r1(indexes), r2(indexes), '*', 'MarkerSize', 6); % Comparison with the Traditional Correlation matrix
% title('Scatter Diagram for Traditional and Weighted Correlations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Weighted Correlation coefficients (uniform weights!!!)', ... % x label ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% ylabel('Traditional Correlation coefficients', ... % y label ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % PLOT HISTOGRAM OF DIFFERENCES BETWEEN COEFFICIENTS
% % OF THE TWO CORRELATION MATRICES
% figure('Position', [50 50 950 630]);
% hist(r2(indexes) - r1(indexes), 100); % Difference between the two correlation matrices
% temp = 'Differences between Traditional and Weighted ';
% temp(2, :) = 'Correlation coefficients (uniform weights!!!)';
% title(temp, 'FontSize', 20, 'FontWeight', 'Bold');
% ylabel('Absolute Frequencies', 'FontSize', 16, 'FontWeight', 'Bold');
% xlabel(temp, 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
%
% % ======================================================================
% % Example 5: more reliable dynamic correlations (it may take some mins):
% % sensitive with respect to recent shocks, not to remote ones!
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 500; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
%
% % CHOOSE WEIGHTS
% delta = 100; % Running window
% alpha = 2 / delta;
% w0 = 1 / sum(exp(((1:delta) - delta) * alpha));
% w = w0 * exp(((1:delta) - delta) * alpha); % weights: exponential decay
%
% % FIND RELEVANT INDEXES AND PLOT A SCATTER DIAGRAM
% I = ones(N);
% indexes = find(tril(I, -1)); % indexes of lower triangular square matrix
%
% % COMPUTE DYNAMIC CORRELATION MATRICES
% for i = 1:(T - delta + 1)
% temp = weightedcorrs(Y(i:(delta + i - 1), :), w); % Dynamic Weighted Correlation matrix
% r1(i) = mean(temp(indexes)); % consider only the lower matrix, in vectorial form, compute the mean
% temp = corrcoef(Y(i:(delta + i - 1), :)); % Dynamic Traditional Correlation matrix
% r2(i) = mean(temp(indexes)); % consider only the lower matrix, in vectorial form, compute the mean
% end
%
% % PLOT THE AVERAGE CORRELATIONS, BOTH WEIGHTED AND TRADITIONAL
% figure('Position', [50 50 950 630]);
% plot(r2, '.g', 'MarkerSize', 24); % the red curve is less sensitive with respect
% hold on; % to remote shocks at time (t - delta + 1) and more
% plot(r1, '.r', 'MarkerSize', 18); % sensitive with respect to present shocks at time t
% title('Moving Average Correlations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Time', 'FontSize', 16, 'FontWeight', 'Bold');
% temp = ' Moving Average Correlations ';
% temp(2, :) = '(green --> traditional; red --> weighted)';
% ylabel(temp, 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
%
% % ======================================================================
% % Example 6: Bell-shaped weights centered on single observations
% % ======================================================================
%
% % GENERATE CORRELATED STOCHASTIC PROCESSES
% N = 2; % number of variables
% T = 1000; % number of observations
% Y = randn(T, N); % shocks from standardized normal distribution
% lambda = 2;
% Y = Y + repmat(poissrnd(lambda, T, 1) .* (2 * rand(T, 1) - 1), 1, N); % shocks plus common factor
% Y = cumsum(Y); % correlated stochastic processes
% alpha = 16 / T / T;
%
% % PLOT BELL-SHAPED WEIGHTS IN THREE CASES: AT THE BEGINNING,
% % IN THE MIDDLE AND AT THE END OF THE PERIOD
% figure('Position', [50 50 950 630]);
% i = 1;
% w = exp(-alpha * ((1:T) - i) .^ 2);
% w = w / sum(w);
% plot(w, 'g', 'LineWidth', 7);
% hold on;
% i = floor(T / 2);
% w = exp(-alpha * ((1:T) - i) .^ 2);
% w = w / sum(w);
% plot(w, 'r', 'LineWidth', 7);
% i = T;
% w = exp(-alpha * ((1:T) - i) .^ 2);
% w = w / sum(w);
% plot(w, 'b', 'LineWidth', 7);
% title('Weights assigned to observations', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Observations', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel('Weights', 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
% legend('Weights centered on t = 1', ... % Legend: first item ...
% sprintf('Weights centered on t = %d', floor(T / 2)), ... % ... second item ...
% sprintf('Weights centered on t = %d', T), ... % ... third item ...
% 'Location', 'North'); % Legend Location
%
% % COMPUTE CORRELATIONS FOR EACH POSSIBLE CENTRAL OBSERVATION
% for i = 1:T
% w = exp(-alpha * ((1:T) - i) .^ 2); % weights: exponential decay
% w = w / sum(w);
% temp = weightedcorrs(Y, w); % Weighted Correlation matrix
% r(i) = temp(2);
% end
% temp = corrcoef(Y);
% rr = temp(2);
%
% % PLOT CENTERED CORRELATION FOR EACH POSSIBLE CENTRAL OBSERVATION
% figure('Position', [50 50 950 630]);
% plot([1, T], [rr, rr], 'b', 'LineWidth', 10);
% hold on;
% plot(r, '.r', 'MarkerSize', 24);
% title('Correlations with bell-shaped weights', ... % title label for the plot ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% xlabel('Central Observation', 'FontSize', 16, 'FontWeight', 'Bold');
% ylabel([' Correlations with bell-shaped weights '; ... % y label, first line ...
% 'centered on different Central Observations'], ... % y label, second line ...
% 'FontSize', 16, 'FontWeight', 'Bold'); % ... font properties
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
% % PLOT HISTOGRAM OF CENTERED CORRELATIONS COMPUTED
% % FOR EACH POSSIBLE CENTRAL OBSERVATION
% figure('Position', [50 50 950 630]);
% hist(r, T / 20);
% title([' Histogram of the correlation computed with '; ... % title label, first line ...
% 'bell-shaped weights on different Central Observations'], ... % title label, second line ...
% 'FontSize', 20, 'FontWeight', 'Bold'); % ... font properties
% ylabel('Absolute Frequencies', 'FontSize', 16, 'FontWeight', 'Bold');
% xlabel('Correlation', 'FontSize', 16, 'FontWeight', 'Bold');
% set(gca, 'FontSize', 14, 'FontWeight', 'Bold');
% grid on;
%
%
% % ======================================================================
%
% See also CORRCOEF, COV, STD, MEAN.
%
% % ======================================================================
%
% Author: Francesco Pozzi
% E-mail: francesco.pozzi@anu.edu.au
% Date: 23 July 2008
%
% % ======================================================================
%
% Check input
ctrl = isvector(w) & isreal(w) & ~any(isnan(w)) & ~any(isinf(w)) & all(w > 0);
if ctrl
w = w(:) / sum(w); % w is column vector
else
error('Check w: it needs be a vector of real numbers with no infinite or nan values!')
end
ctrl = isreal(Y) & ~any(isnan(Y)) & ~any(isinf(Y)) & (size(size(Y), 2) == 2);
if ~ctrl
error('Check Y: it needs be a 2D matrix of real numbers with no infinite or nan values!')
end
ctrl = length(w) == size(Y, 1);
if ~ctrl
error('size(Y, 1) has to be equal to length(w)!')
end
[T, N] = size(Y); % T: number of observations; N: number of variables
temp = Y - repmat(w' * Y, T, 1); % Remove mean (which is, also, weighted)
temp = temp' * (temp .* repmat(w, 1, N)); % Covariance Matrix (which is weighted)
temp = 0.5 * (temp + temp'); % Must be exactly symmetric
R = diag(temp); % Variances
R = temp ./ sqrt(R * R'); % Matrix of Weighted Correlation Coefficients
