-
Notifications
You must be signed in to change notification settings - Fork 2
/
RunMeExampleLocallyAdaptiveEDS.m
123 lines (102 loc) · 5.82 KB
/
RunMeExampleLocallyAdaptiveEDS.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
% RunMeExampleLocallyAdaptiveEDS.m is a simple example that illustrates
% the construction of locally adaptive epsilon-distinguishable sets, as
% described in the article "Merging simulation and projection approaches
% to solve high-dimensional problems with an application to a new Keynesian
% model" by Lilia Maliar and Serguei Maliar (2015), Quantitative Economics
% 6/1, pages 1–47 (henceforth, MM, 2015).
%
% This version: March 19, 2015. First version: June 27, 2011.
%
% -------------------------------------------------------------------------
% Copyright © 2015 by Lilia Maliar and Serguei Maliar. All rights reserved.
% The code may be used, modified and redistributed under the terms provided
% in the file "License_Agreement.txt".
% -------------------------------------------------------------------------
clc;
clear all;
% 1. Parameters of the model
% ---------------------------
alpha = 0.36; % Capital share
beta = 0.99; % Discount factor
d = 1; % Depreciation rate
sigma = 0.01; % Standard deviation for log shock
rho = 0.95; % Persistence of log technology shock
A = (1/beta-(1-d))/alpha; % Technology level
% 2. Simulate time-series solution to the neoclassical growth model
% -----------------------------------------------------------------
% Under the assumption of full depreciation of capital and logarithmic utility
% function, this model admits a closed form solution
% 2.1 Allocate memory for the simulated series
% --------------------------------------------
T = 1000000; % Length of simulation
a = ones(T,1); % Productivity levels
k = ones(T,1); % Capital
% 2.2 Simulate capital and productivity series
% --------------------------------------------
epsi = randn(T,1)*sigma; % Draw productivity shocks
for t = 1:T-1;
k(t+1,1) = alpha*beta*A*a(t)*k(t)^alpha; % Capital
a(t+1,:) = a(t,:).^rho.*exp(epsi(t,:)); % Productivity
end;
% 3. Reduce the number of simulated points by selecting each 100th point
% -----------------------------------------------------------------------
% This is a cheap method that allows to obtain simulated points that are
% uncorrelated (independent) across time
N = 10000; % Select N = 10,000 points out of T = 1,000,0000 points
a1 = ones(N,1); % Allocate memory for productivity
k1 = ones(N,1); % Allocate memory for capital
for j = 1:10000
a1(j,:) = a(j*100,1);
% Select each 100th point
k1(j,:) = k(j*100,1);
end
Data = [k1 a1]; % Data that will be used to construct the EDS grid
% 4. Estimating density function on principal components (PCs)
% ------------------------------------------------------------
[density PCn] = Density(Data,Data); % Estimate the density function in
% each data point
Sorted_PCn=sortrows([PCn density],3); % Sort principal components by
% density function
Sorted_Data=sortrows([Data density],3); % Sort Data by density function
cutoff = 500; % The number of lowest density
% points to be dropped
Data1 = Sorted_Data(1+cutoff:end,1:2); % Data after removing the lowest
% density points
PCn1 = Sorted_PCn(1+cutoff:end,1:2); % Principal components after
% removing the lowest density
% points
n = size(PCn1,1); % Number of points in PCn1
% 5. Construct locally adaptive EDS
% ---------------------------------
vepsilon = 0.3*ones(n,1).*(PCn1(:,1)>0)+1.2*ones(n,1).*(PCn1(:,1)<=0);
% define vector of epsilon of
% length of PC1 and set epsilon=0.5
% if PCn1>0, and set epsilon=1.2
% if PCn1<=0; this yields more
% dense grid when PCn1>0 than when
% PCn1<=0
[EDS,EDS_PCn] = LocallyAdaptiveEDS(Data1,vepsilon);
% Construct EDS set
N_EDS = size(EDS,1); % Exact number of points in the EDS set
% 6. Figures
% ----------
subplot(3,2,1);
plot(k1(:,1),a1(:,1),'.','MarkerSize',2), xlabel('k_t'), ylabel('a_t'),axis([0.82 1.18 0.85 1.15]);
title('Figure 1a. Simulated points');
subplot(3,2,2);
plot(Sorted_PCn(:,1),Sorted_PCn(:,2),'.','MarkerSize',2), xlabel('PC^1_t'), ylabel('PC^2_t'),axis([-5.5 5.5 -5.5 5.5]);
title('Figure 1b. Principal components (PCs)');
subplot(3,2,3);
axis([-4 4 -4 4]);
plot(Sorted_PCn(cutoff+1:10000,1),Sorted_PCn(cutoff+1:10000,2),'.',Sorted_PCn(1:cutoff,1),Sorted_PCn(1:cutoff,2),'x','MarkerSize',4), xlabel('PC^1_t'), ylabel('PC^2_t'),axis([-5.5 5.5 -5.5 5.5]);
title('Figure 1c. Dropped off low density points');
subplot(3,2,4);
plot(Sorted_PCn(cutoff+1:10000,1),Sorted_PCn(cutoff+1:10000,2),'.',EDS_PCn(:,1),EDS_PCn(:,2),'ro','MarkerSize',4'), xlabel('PC^1_t'), ylabel('PC^2_t'),axis([-5.5 5.5 -5.5 5.5]);
title('Figure 1d. EDS balls on PCs');
subplot(3,2,5);
plot(EDS_PCn(:,1),EDS_PCn(:,2),'o','MarkerSize',4'), xlabel('PC^1_t'), ylabel('PC^2_t'),axis([-5.5 5.5 -5.5 5.5]);
title('Figure 1e. EDS on PCs');
subplot(3,2,6);
axis([-4 4 -4 4]);
plot(EDS(:,1),EDS(:,2),'o','MarkerSize',4), xlabel('k_t'), ylabel('a_t'),axis([0.82 1.18 0.85 1.15]);
title('Figure 1f. EDS on original data');