% DELTA_AS.M            05-31-2011         code by: Komunjer and Ng
% This function computes the delta matrix in Proposition 2-S (SINGULAR CASE)
% of Komunjer and Ng (2011). The derivatives are computed numerically.

% Adjusted by Filippo Ferroni 05-31-2012


function [Delta,lambda_grad,T_grad,U_grad] = delta_as(theta,A,B,C,D,Sigma,variabs)

n_x    = size(A,1);
n_eps  = size(B,2);
n_y    = size(C,1);
lambda = [vec(A); vec(B); vec(C); vec(D); vec(Sigma)];

% computes Delta_theta
lambda_grad = zeros(size(lambda,1),size(theta,1));
for i=1:1:size(theta)
    delta_theta = zeros(size(theta)); 
%     delta_theta(i) = theta(i)*1e-3; % original KN(2011)
    delta_theta(i) = theta(i)*1e-2;
    if delta_theta(i)==0; delta_theta(i) = 1e-2; end;
    theta_p = theta + delta_theta; 
    [minA,minB,minC,minD,Sigma]=sw4_abcd(theta_p,variabs);
%     [minA,minB,minC,minD]=put2minimal(n_x,A_0,B_0,C_0);
    lambda_p = [vec(minA); vec(minB); vec(minC); vec(minD); vec(Sigma)];
    theta_m = theta - delta_theta; 
    [minA,minB,minC,minD,Sigma]=sw4_abcd(theta_m,variabs);
%    [minA,minB,minC,minD]=put2minimal(n_x,A_0,B_0,C_0);
    lambda_m = [vec(minA); vec(minB); vec(minC); vec(minD); vec(Sigma)];
    lambda_grad(:,i) = (lambda_p - lambda_m)/(2*delta_theta(i));
end;

% computes the permutation matix T
T = [];
for j=1:1:n_eps
    ind_j = zeros(n_eps,1); ind_j(j) = 1;
    T = [T, kron(eye(n_eps,n_eps),ind_j)];
end;

% computes Delta_T
T_grad = [kron(A',eye(n_x)) - kron(eye(n_x),A);
          kron(B',eye(n_x));
          -1*kron(eye(n_x),C);
          zeros(n_y*n_eps,n_x^2);
          zeros(n_eps^2,n_x^2)];

%computes Delta_U
U_grad = [zeros(n_x^2,n_eps^2); 
    kron(eye(n_eps),B); 
    zeros(n_y*n_x,n_eps^2); 
    kron(eye(n_eps),D);
    -1*(eye(n_eps^2) + T)*kron(Sigma,eye(n_eps))];

 Delta = [lambda_grad  T_grad   U_grad ];

