function [w,beta,theta,rho,wvec,betavec,thetavec,rhovec] = GenerateMatrices(n)
% generate "reasonable" model parameters w_{ij}, \beta_{ij}, theta_{ij}, \rho_{ij}
% for a large bank like JP Morgan for N=17 sectors and 
% specified n (# companies per sector). The w_{ij} are based on call report
% data

% N     Loan type/sector                            
% 1     construction and land development               
% 2     farmland                                        
% 3     1-4 residential properties                      
% 4     5+ residential properties                       
% 5     nonfarm, non-residential                        
% 6     commercial US banks                             
% 7     banks in foreign countries                      
% 8     agricultural loans                              
% 9     US commercial/industrial loans                  
% 10	non-US commercial/industrial loans              
% 11	Credit Cards                                    
% 12	Other revolving credit plans                    
% 13	automobile loans                                
% 14	other consumer loans                            
% 15	foreign governments                             
% 16	US states and subdivisions                      
% 17	non-depository financial institutions & other   

wvec = [      
   0.006038934108343
   0.000334001336005
   0.321727182643328
   0.071374977443038
   0.042981697993143
   0.004997356482317
   0.029969892391418
   0.001247360439679
   0.143856433245638
   0.053224933278880
   0.041455739282673
   0.004090329157525
   0.065719115957033
   0.031400874418664
   0.001632016006737
   0.020071739054729
   0.159877416760852];

rhovec = [     
    0.3300
    0.3200
    0.3500
    0.3100
    0.3700
    0.9500
    0.9300
    0.5200
    0.4800
    0.5000
    0.1500
    0.1700
    0.2100
    0.1600
    0.9200
    0.9700
    0.8500];
       
% rhovec = [
%     0.23
%     0.22
%     0.25
%     0.21
%     0.27
%     0.15
%     0.13
%     0.12
%     0.18
%     0.28
%     0.15
%     0.17
%     0.21
%     0.16
%     0.12
%     0.17
%     0.15
%     ];

thetavec = [ 
  -2.326347874040841
  -2.326347874040841
  -2.326347874040841
  -2.326347874040841
  -2.326347874040841
  -3.090232306167814
  -3.090232306167814
  -2.326347874040841
  -2.326347874040841
  -2.326347874040841
  -1.644853626951473
  -1.644853626951473
  -1.644853626951473
  -1.644853626951473
  -3.719016485455709
  -3.719016485455709
  -2.326347874040841];
       
betavec = [        
    0.8000
    0.9000
    0.0500
    0.1000
    0.1000
    0.9000
    0.1000
    0.7000
    0.0500
    0.1000
    0.1000
    0.8000
    0.1000
    0.1000
    0.8000
    0.1000
    0.0500];

N=length(wvec);
w = zeros(N,n);
rho = zeros(N,n);
theta = zeros(N,n);
for i=1:N
    w(i,:) = wvec(i)/n;
    rho(i,:) = rhovec(i) + 0.01*randn(1,n);
    theta(i,:) = thetavec(i) + 0.01*randn(1,n);
    beta(i,:) = betavec(i) + 0.01*rand(1,n);
end



function [VaR,VaR_upper,VaR_lower,R_pi] = MCsimulatedVaR(rho,beta,theta,w,c,q,num_trials,confidence,qs)
%
% This function produces a Monte Carlo simulated VaR_q plot for the 
% hierarchical multi-factor credit
% portfolio problem and generates error bars at values of q specified
% in the qs vector
%
% rho:     matrix of asset return variances
% beta:    matrix of factor loadings
% theta:   matrix of threshold default values s.t. prob default =
%          Phi(theta_{ij})
% w:       matrix of exposures
% c:       loss given default (LGD). It's a scalar e.g. c=1.
% q:       vector of risk-values at which to evaluate VaR. Must be between
%          0 and 1 e.g. q = linspace(0,1,50)
%
% num_trials:            number of portfolios to generate e.g. num_trials = 5000
% confidence:            confidence values to create error bars. It's a
%                        number between 0 (no confidence) and 1 (complete confidence).
% VaR:                   Values-at-Risk corresponding to q
% [VaR_upper,VaR_lower]: The upper and lower error bars for VaR,
%                        calculated using the bootstrap method
% qs:                    Values of q at which to calculate error bars e.g.
%                        qs = [0.2 0.4 0.6 0.8]
% R_pi:                  Portfolio loss for each of the num_trials trials

NNN = 50000; % number of bootstrap samples
beta_hat = sqrt(1-beta.^2);

% parameter values
[N,n] = size(rho); % N: number of sectors, n: number of companies per sector
R_pi = zeros(1,num_trials);
for i=1:num_trials
    R_pi(i) = one_draw(rho,beta,beta_hat,theta,w,c,N,n);
end

P = zeros(NNN,length(qs));
for i=1:NNN
    if mod(i,2000) == 0
        sprintf('Generating bootstrap samples: %d/%d',i,NNN)
    end
    % generate num_trials random integers from [1:num_trials] 
    k = randi(num_trials,1,num_trials);
    bootstrap_sample = R_pi(k);
    P(i,:) = quantile(bootstrap_sample,qs);
end
alpha = (1-confidence)/2; % e.g. alpha = 0.05 for confidence_percentage = 0.9
VaR_lower = quantile(P,alpha);
VaR_upper = quantile(P,1-alpha);
VaR = quantile(R_pi,q); VaR(end) = 1;
semilogy(q,VaR,'r--','LineWidth',2);
hold on; 

for i=1:length(qs)
    val(i) = (VaR_lower(i) + VaR_upper(i))/2;
    errorbar(qs(i),val(i),val(i)-VaR_lower(i),VaR_upper(i)-val(i),'r-','LineWidth',2);
end

axis([0 1 1e-3 1e-1]);
h=gca;
set(h,'FontSize',14,'FontName','Times');
end

function R_pi = one_draw(rho,beta,beta_hat,theta,w,c,N,n)
    % different rows of rho, beta etc. correspond to different sectors, 
    % following the convention in the paper
    epsilon_hat = randn(1,1);
    epsilon = kron(randn(N,1),ones(1,n));
    Z1 = sqrt(rho) .* beta_hat*epsilon_hat;
    Z2 = sqrt(rho) .* beta .* epsilon;
    Z3 = sqrt(1-rho) .* randn(N,n);
    Z = Z1+Z2+Z3;
    LOSSES = c*(Z <= theta);
    R_pi = sum(sum( w.*LOSSES ));
end

function VaR = AnalyticVaR(rho,beta,theta,w,c,q)
% provides analytic solution for N sectors and n companies per sector
% rho = A(:,1); beta = A(:,2); theta = A(:,3); w = A(:,4);
beta_hat = sqrt(1-beta.^2);
[N,n] = size(rho);
for i=1:length(q)
    arg = ( theta + sqrt(rho).*beta_hat.*invPhi(q(i)) )./ (sqrt(1-rho+beta.^2.*rho));
    VaR(i) = c*sum(sum(w.* Phi(arg)));
end

semilogy(q,VaR,'b-','LineWidth',2); hold on;
xlabel('q','FontSize',14,'FontName','Times','FontAngle','italic');
ylabel('VaR_q','FontSize',14,'FontName','Times');
end


function out = Phi(x)
out = 1/2*(1+erf(x/sqrt(2)));
end

function out = invPhi(x)
out = sqrt(2)*erfinv(2*x-1);
end