function [X, val_obj] = Correlation_Newton(G,C,m,n,tau,tau0,tol)
%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% This code is for the Local Correlation Stress Testing of the Lagrangian Dual approach described in
%%%%%% Houduo Qi and Defeng Sun, "Correlation Stress Testing for Value-at-Risk: An Unconstrained
%%%%%% Convex Optimization Approach", Department of Mathematics, National
%%%%%% University of Singapore, March 2007
%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%
% for solving
%% min 0.5 || X - C ||^2
%% s.t. X_ii = 1, i=1, ..., n
%% X_ij = C_ij, i=1, ..., m; j =1, ..., n (Recall X = X^T)
%% X- tau*I >=0 (symmetric and positive semi-definite)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%
% G: the unstresed correlation matrix (n by n)
% C: stressed correlation matrix (n by n)
% C(1:m,1:m) is fixed
% tau0: the threshold for positive definiteness of C
% tau: the lower bound for the smallest eigevalue of X
% tol: stopping crierion for the KKT system: 1.0e-5 ~ 1.0e-6
% X: the calibrated correlation matrix
% val_obj: final objective function value
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Last updated on March 20, 2007 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ---Newton method starts--- ')
t0=cputime;
G =(G+G')/2; % make G symmetric
C =(C+C')/2;
if m<1
tau =tau0;
else
tau =min(tau, 0.90*tau0);
end
n_2=n-m;
if m>0
rhs1 =ones(n,1);
[P,D] = eig(G);
d = diag(D);
d =real(d);
if min(d) < tau0
fprintf('\n')
fprintf('The smallest eignevalue of the initial G < %d \n', tau0)
disp(' ---Start the Pre-processing--- ')
[G,y]=CorNewton1(G,rhs1,tau0);
C(1:m,1:n) = G(1:m,1:n); % update G to make sure G is positive definite
C((m+1):m,1:m) = G((m+1):m,1:m);
disp(' ---Pre-processing finished--- ')
else
fprintf('\n')
fprintf('The smallest eignevalue of the initial matrix G(1:m,1:m) >= %d \n', tau0)
disp(' ---No pre-processing needed--- ')
end
end
C =C-tau*eye(n); % reset C: C>= tau *I;
C_1=zeros(m, m);
C_2=zeros(m, n_2);
C_1=C(1:m, 1:m);
C_2=C(1:m, (m+1):n);
% Calculate Schur Complement
Y=C_2'*inv(C_1)*C_2;
dd=diag(Y);
C_3=zeros(n_2, n_2);
C_3=C((m+1):n, (m+1):n);
b=(1-tau)*ones(n_2,1)-dd;
C_3 = C_3-Y; % C_3 updated
[Y0,y0] = CorNewton1(C_3,b,0); % updated C_3 is calibrated
C_3 = Y0+Y; % C_3 is updated back
X = [C_1,C_2; C_2',C_3];
G = X-C;
val_obj = sum(sum(G.*G))/2;
X = X + tau*eye(n); % optimal X
fprintf('Newton: Final function value %d \n', val_obj)
time_used = cputime -t0