-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathfs_cholesky_decomposition_d5.m
153 lines (144 loc) · 4.89 KB
/
fs_cholesky_decomposition_d5.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
%% Cholesky Decomposition Correlated Five Dimensional Multivariate Normal Shock
% *back to* <https://fanwangecon.github.io *Fan*>*'s* <https://fanwangecon.github.io/Math4Econ/
% *Intro Math for Econ*>*,* <https://fanwangecon.github.io/M4Econ/ *Matlab Examples*>*,
% or* <https://fanwangecon.github.io/MEconTools/ *MEconTools*> *Repositories*
%
% Generate variance-covariance matrix from correlation and standard deviation.
% Draw five correlated normal shocks using the MVNRND function. Draw five correlated
% normal shocks from uniform random variables using Cholesky Decomposition.
%%
% * <https://fanwangecon.github.io/M4Econ/simulation/normal/htmlpdfm/fs_cholesky_decomposition.html
% fs_cholesky_decomposition>
% * <https://fanwangecon.github.io/M4Econ/simulation/normal/htmlpdfm/fs_cholesky_decomposition_d5.html
% fs_cholesky_decomposition_d5>
% * <https://fanwangecon.github.io/M4Econ/simulation/normal/htmlpdfm/fs_bivariate_normal.html
% fs_bivariate_normal>
%% Correlation and Standard Deviations to Variance Covariance Matrix
% Given correlations and standard deviations, what is the variance and covariance
% matrix? Assume mean 0. The first three variables are correlated, the final two
% are iid.
%
% First the ingredients:
% mean array
ar_mu = [0,0,0,0,0];
% standard deviations
ar_sd = [0.3301, 0.3329, 0.3308, 2312, 13394];
% correlations
mt_cor = ...
[1,0.1226,0.0182,0,0;...
0.1226,1,0.4727,0,0;...
0.0182,0.4727,1,0,0;...
0,0,0,1,0;...
0,0,0,0,1];
% show
disp(mt_cor);
%%
% Second, we know that variance is the square of standard deviation. And we
% know the formula for covariance, which is variance divided by two standard deviations.
% So for the example here, we have:
% initialize
mt_varcov = zeros([5,5]);
% variance
mt_varcov(eye(5)==1) = ar_sd.^2;
% covariance
mt_varcov(1,2) = mt_cor(1,2)*ar_sd(1)*ar_sd(2);
mt_varcov(2,1) = mt_varcov(1,2);
mt_varcov(1,3) = mt_cor(1,3)*ar_sd(1)*ar_sd(3);
mt_varcov(3,1) = mt_varcov(1,3);
mt_varcov(2,3) = mt_cor(2,3)*ar_sd(2)*ar_sd(3);
mt_varcov(3,2) = mt_varcov(2,3);
% show
disp(mt_varcov(1:3,1:3));
disp(mt_varcov(4:5,4:5));
%% Draw Five Correlated Shocks Using MVNRND
% Generate N5 correlated shock structure
% Generate Scores
rng(123);
N = 50000;
mt_kw97_eps = mvnrnd(ar_mu, mt_varcov, N);
% graph
figure();
% subfigure 1
subplot(2,2,1);
scatter(mt_kw97_eps(:,1), mt_kw97_eps(:,2));
ylabel('White Collar Wage Shock');
xlabel('Blue Collar Wage Shock')
grid on;
% subfigure 2
subplot(2,2,2);
scatter(mt_kw97_eps(:,1), mt_kw97_eps(:,3));
ylabel('White Collar Wage Shock');
xlabel('Military Wage Shock')
grid on;
% subfigure 3
subplot(2,2,3);
scatter(mt_kw97_eps(:,3), mt_kw97_eps(:,2));
ylabel('Military Wage Shock');
xlabel('Blue Collar Wage Shock')
grid on;
% subfigure 4
subplot(2,2,4);
scatter(mt_kw97_eps(:,1), mt_kw97_eps(:,4));
ylabel('White Collar Wage Shock');
xlabel('School Shock')
grid on;
%%
% What are the covariance and correlation statistics?
disp([num2str(round(corrcoef(mt_kw97_eps),3))]);
disp([num2str(round(corrcoef(mt_kw97_eps),2))]);
%% Draw Five Correlated Shocks Using Cholesky Decomposition
% Following what we did for bivariate normal distribution, we can now do the
% same for five different shocks at the same time (For more details see <https://eml.berkeley.edu/~train/distant.html
% Train (2009)>):
%%
% # Draw 5 normal random variables that are uncorrelated
% # Generate 5 correlated shocks
%%
% Draw the shocks
% Draws uniform and invert to standard normal draws
rng(123);
ar_unif_draws = rand(1,N*5);
ar_normal_draws = norminv(ar_unif_draws);
ar_draws_eta_1 = ar_normal_draws((N*0+1):N*1);
ar_draws_eta_2 = ar_normal_draws((N*1+1):N*2);
ar_draws_eta_3 = ar_normal_draws((N*2+1):N*3);
ar_draws_eta_4 = ar_normal_draws((N*3+1):N*4);
ar_draws_eta_5 = ar_normal_draws((N*4+1):N*5);
% Choesley decompose the variance covariance matrix
mt_varcov_chol = chol(mt_varcov, 'lower');
% Generate correlated random normals
mt_kp97_eps_chol = ar_mu' + mt_varcov_chol*([...
ar_draws_eta_1; ar_draws_eta_2; ar_draws_eta_3; ar_draws_eta_4; ar_draws_eta_5]);
mt_kp97_eps_chol = mt_kp97_eps_chol';
%%
% Graph:
% graph
figure();
% subfigure 1
subplot(2,2,1);
scatter(mt_kp97_eps_chol(:,1), mt_kp97_eps_chol(:,2));
ylabel('CHOL White Collar Wage Shock');
xlabel('CHOL Blue Collar Wage Shock')
grid on;
% subfigure 2
subplot(2,2,2);
scatter(mt_kp97_eps_chol(:,1), mt_kp97_eps_chol(:,3));
ylabel('CHOL White Collar Wage Shock');
xlabel('CHOL Military Wage Shock')
grid on;
% subfigure 3
subplot(2,2,3);
scatter(mt_kp97_eps_chol(:,3), mt_kp97_eps_chol(:,2));
ylabel('CHOL Military Wage Shock');
xlabel('CHOL Blue Collar Wage Shock')
grid on;
% subfigure 4
subplot(2,2,4);
scatter(mt_kp97_eps_chol(:,1), mt_kp97_eps_chol(:,4));
ylabel('CHOL White Collar Wage Shock');
xlabel('CHOL School Shock')
grid on;
%%
% What are the covariance and correlation statistics?
disp([num2str(round(corrcoef(mt_kp97_eps_chol),3))]);
disp([num2str(round(corrcoef(mt_kp97_eps_chol),2))]);