Add files via upload

rtol · Apr 4, 2023 · 6292675 · 6292675
1 parent 4d662ce
commit 6292675
Show file tree

Hide file tree

Showing 14 changed files with 983 additions and 237 deletions.
diff --git a/BVDecompTest.m b/BVDecompTest.m
@@ -0,0 +1,39 @@
+function [chi, p] = BVDecompTest(K,NQ,N)
+%function [chi, p] = DecompositionTest(K,NQ,N)
+%tests whether the components of a kernel deviates from the composite
+%
+%K is the matrix of component kernels, from
+%[K,grid] = KernelDecomposition(data,partition,kernel)
+%NQ>1 is the number of quantiles
+%N is the number of data used to estimate the kernels
+%
+%chi is the test-statistic, distributed as a chi-squared with degrees of
+%freedom equal to the number of quantiles minus one times the number of
+%component minus one
+%p is the p-value
+%
+%First version, Richard Tol, 11 March 2020
+%This version, Richard Tol, 11 March 2020
+P = size(K,2);
+G = size(K,1);
+CDF = K(:,1);
+for i=2:G,
+    CDF(i) = CDF(i-1)+K(i,1);
+end
+
+for i=1:NQ,
+    Q(i,:) = sum(K(CDF<=i/NQ & CDF > (i-1)/NQ,2:P));
+end
+
+w = sum(K(:,2:P),1);
+
+for i=1:NQ,
+    Pearson(i,:) = ((Q(i,:) - w/NQ).^2)./(w/NQ);
+end
+
+chi = sum(sum(Pearson))*N/NQ;
+
+p = 1-chi2cdf(chi,(P-2)*(NQ-1));
+
+end
+
diff --git a/Characteristics.m b/Characteristics.m
@@ -1,50 +1,55 @@
-%% Characteristics
-%
-% First version: Richard Tol, 2 November 2011
-% This version: Richard Tol, 2 November 2011
-
-display('Characteristics');
-
-for i=1:NFilter,
-    KernelMean(i) = SCCgrid*JointPDF(:,i);
-    vmeansq= (SCCgrid.*SCCgrid)*JointPDF(:,i);
-    var = vmeansq-KernelMean(i)*KernelMean(i);
-    KernelStDev(i) = sqrt(var);
-    [vmax vind] = max(JointPDF(:,i));
-    KernelMode(i) = SCCgrid(vind);
-    vsmaller = (JointCDF(:,i) < 0.5);
-    vlarger = (JointCDF(:,i) > 0.5);
-    KernelMedian(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-    vsmaller = (JointCDF(:,i) < 0.05);
-    vlarger = (JointCDF(:,i) > 0.05);
-    Kernel05(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-    vsmaller = (JointCDF(:,i) < 0.95);
-    vlarger = (JointCDF(:,i) > 0.95);
-    Kernel95(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-end
-
-
-%{
-for i=1:NTime,
-    for f=1:NFilter,
-        if sum(CondCDF(:,i,f)) > 0,
-            TimeMean(i,f) = SCCgrid*CondPDF(:,i,f);
-            vmeansq = (SCCgrid.*SCCgrid)*CondPDF(:,i,f);
-            var = vmeansq - TimeMean(i,f)*TimeMean(i,f);
-            TimeStDev(i,f) = sqrt(var);
-            [vmax vind] = max(CondPDF(:,i,f));
-            TimeMode(i,f) = SCCgrid(vind);
-            vsmaller = (CondCDF(:,i,f) < 0.5);
-            vlarger = (CondCDF(:,i,f) > 0.5);
-            TimeMedian(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-            vsmaller = (CondCDF(:,i,f) < 0.05);
-            vlarger = (CondCDF(:,i,f) > 0.05);
-            Time05(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-            vsmaller = (CondCDF(:,i,f) < 0.95);
-            vlarger = (CondCDF(:,i,f) > 0.95);
-            Time95(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
-        end
-    end
-end
-%}    
+%% Characteristics
+%
+% First version: Richard Tol, 2 November 2011
+% This version: Richard Tol, 2 November 2011
+
+display('Characteristics');
+
+for i=1:NFilter,
+    KernelMean(i) = SCCgrid*JointPDF(:,i);
+    vmeansq= (SCCgrid.*SCCgrid)*JointPDF(:,i);
+    var = vmeansq-KernelMean(i)*KernelMean(i);
+    KernelStDev(i) = sqrt(var);
+    [vmax vind] = max(JointPDF(:,i));
+    KernelMode(i) = SCCgrid(vind);
+    vsmaller = (JointCDF(:,i) < 0.5);
+    vlarger = (JointCDF(:,i) > 0.5);
+    KernelMedian(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+    vsmaller = (JointCDF(:,i) < 0.05);
+    vlarger = (JointCDF(:,i) > 0.05);
+    Kernel05(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+    vsmaller = (JointCDF(:,i) < 0.95);
+    vlarger = (JointCDF(:,i) > 0.95);
+    Kernel95(i) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+end
+
+%% kernelmeans.tex
+Report = [KernelMean(18:23);KernelMean(24:29);KernelMean(30:35);KernelMean(36:41);KernelMean(42:47)];
+Report2 = [KernelMean(1) KernelMean(2) KernelMean(3) KernelMean(5) KernelMean(7);
+ KernelStDev(1) KernelStDev(2) KernelStDev(3) KernelStDev(5) KernelStDev(7);
+ KernelMedian(1) KernelMedian(2) KernelMedian(3) KernelMedian(5) KernelMedian(7);
+ KernelMode(1) KernelMode(2) KernelMode(3) KernelMode(5) KernelMode(7)];
+%{
+for i=1:NTime,
+    for f=1:NFilter,
+        if sum(CondCDF(:,i,f)) > 0,
+            TimeMean(i,f) = SCCgrid*CondPDF(:,i,f);
+            vmeansq = (SCCgrid.*SCCgrid)*CondPDF(:,i,f);
+            var = vmeansq - TimeMean(i,f)*TimeMean(i,f);
+            TimeStDev(i,f) = sqrt(var);
+            [vmax vind] = max(CondPDF(:,i,f));
+            TimeMode(i,f) = SCCgrid(vind);
+            vsmaller = (CondCDF(:,i,f) < 0.5);
+            vlarger = (CondCDF(:,i,f) > 0.5);
+            TimeMedian(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+            vsmaller = (CondCDF(:,i,f) < 0.05);
+            vlarger = (CondCDF(:,i,f) > 0.05);
+            Time05(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+            vsmaller = (CondCDF(:,i,f) < 0.95);
+            vlarger = (CondCDF(:,i,f) > 0.95);
+            Time95(i,f) = (min(SCCgrid(vlarger)) + max(SCCgrid(vsmaller)))/2;
+        end
+    end
+end
+%}    
 clear v*
diff --git a/ConfirmationBias.m b/ConfirmationBias.m
@@ -1,44 +1,44 @@
-%% ConfirmationBias - Tests for confirmation bias
-%
-% First version: Richard Tol, 1 November 2011
-% This version: Richard Tol, 4 April 2015
-
-display('Test confirmation bias');
-
-StartYear = 1991;
-EndYear = Newest;
-BiasYear(1)= StartYear;
-
-for j=1:NFilter,
-    clear Average StDev Max Min Mean;
-    FilterYear = (Year <= 1991) & Filter(:,j);
-    Average(1) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
-    AveSq = sum(SCC(FilterYear).*SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
-    StDev(1) = sqrt(AveSq-Average(1)*Average(1));
-
-    for i= 2:(EndYear - StartYear + 1),
-        BiasYear(i) = StartYear+i-1;
-        FilterYear = (Year <= BiasYear(i)) & Filter(:,j);
-        Average(i) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
-        AveSq = sum(SCC(FilterYear).*SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
-        StDev(i) = sqrt(AveSq-Average(i)*Average(i));
-        FilterYear = (Year == BiasYear(i)) & Filter(:,j);
-        Mean(i) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
-        if length(SCC(FilterYear)) > 0,
-            Max(i) = max(SCC(FilterYear));
-            Min(i) = min(SCC(FilterYear));
-        else
-            Max(i) = NaN;
-            Min(i) = NaN;
-        end
-    end
-
-    subplot(NFilter,1,j);
-    plot(BiasYear,Average,'b',BiasYear,Average-2*StDev,'b:',BiasYear,Average+2*StDev,'b:',BiasYear,Mean,'r*',BiasYear,Min,'rx',BiasYear,Max,'rx');
-    title(Titles(j,:));
-    xlabel('year of publication')
-    ylabel('dollar per metric ton of carbon');
-    legend('average of previous studies','upper bound of previous studies','lower bound of previous studies','mean of current studies','maximum of current studies','minimum of current studies');
-end
-
+%% ConfirmationBias - Tests for confirmation bias
+%
+% First version: Richard Tol, 1 November 2011
+% This version: Richard Tol, 4 April 2015
+
+display('Test confirmation bias');
+
+StartYear = 1991;
+EndYear = Newest;
+BiasYear(1)= StartYear;
+
+for j=1:NFilter,
+    clear Average StDev Max Min Mean;
+    FilterYear = (Year <= 1991) & Filter(:,j);
+    Average(1) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
+    AveSq = sum(SCC(FilterYear).*SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
+    StDev(1) = sqrt(AveSq-Average(1)*Average(1));
+
+    for i= 2:(EndYear - StartYear + 1),
+        BiasYear(i) = StartYear+i-1;
+        FilterYear = (Year <= BiasYear(i)) & Filter(:,j);
+        Average(i) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
+        AveSq = sum(SCC(FilterYear).*SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
+        StDev(i) = sqrt(AveSq-Average(i)*Average(i));
+        FilterYear = (Year == BiasYear(i)) & Filter(:,j);
+        Mean(i) = sum(SCC(FilterYear).*AuthorWeight(FilterYear))/sum(AuthorWeight(FilterYear));
+        if length(SCC(FilterYear)) > 0,
+            Max(i) = max(SCC(FilterYear));
+            Min(i) = min(SCC(FilterYear));
+        else
+            Max(i) = NaN;
+            Min(i) = NaN;
+        end
+    end
+
+    subplot(NFilter,1,j);
+    plot(BiasYear,Average,'b',BiasYear,Average-2*StDev,'b:',BiasYear,Average+2*StDev,'b:',BiasYear,Mean,'r*',BiasYear,Min,'rx',BiasYear,Max,'rx');
+    title(Titles(j,:));
+    xlabel('year of publication')
+    ylabel('dollar per metric ton of carbon');
+    legend('average of previous studies','upper bound of previous studies','lower bound of previous studies','mean of current studies','maximum of current studies','minimum of current studies');
+end
+
 clear v*
diff --git a/ConstructBivarPDF.m b/ConstructBivarPDF.m
@@ -1,37 +1,37 @@
-%% Construct PDF Bivariate
-%
-% First version: Richard Tol, 30 October 2012
-% This version: Richard Tol, 30 October 2012
-
-display('Construct bivariate PDF');
-
-TimeGrid(1) = Oldest;
-NTime = Newest-Oldest+1;
-for i=2:NTime,
-    TimeGrid(i)=TimeGrid(i-1)+1;
-end
-
-%Bivar = JointPDF(:,1)*PDFTime;
-
-Bivar = zeros(NGrid+1,NTime,NFilter);
-
-for i=1:NGrid+1,
-    for j=1:NTime,
-        for f=1:NFilter,
-            for k=1:NEstimates,
-                %Bivar(i,j) = Bivar(i,j) + exp(-(SCCgrid(i)-SCC(k))/SampleStDev(1)).*exp(-exp(-(SCCgrid(i)-SCC(k))/SampleStDev(1)))/SampleStDev(1) * Weight(k,1) * exp(-0.5*((TimeGrid(j)-Year(k))/StDevYear)^2);
-                if Year(k) >= TimeGrid(j),
-                   Bivar(i,j,f) = Bivar(i,j,f) + exp(-(SCCgrid(i)-SCC(k))/SampleStDev(f)).*exp(-exp(-(SCCgrid(i)-SCC(k))/SampleStDev(f)))/SampleStDev(f) * Weight(k,f) * exp(-(Year(k)-TimeGrid(j)));
-                end
-            end
-        end
-    end
-end
-for f = 1:NFilter,
-    vsum= sum(sum(Bivar(:,:,f)));
-    Bivar(:,:,f) = Bivar(:,:,f)/vsum;
-    MargTime(:,f) = sum(Bivar(:,:,f),1);
-    MargSCC(:,f) = sum(Bivar(:,:,f),2);
-end
-
+%% Construct PDF Bivariate
+%
+% First version: Richard Tol, 30 October 2012
+% This version: Richard Tol, 30 October 2012
+
+display('Construct bivariate PDF');
+
+TimeGrid(1) = Oldest;
+NTime = Newest-Oldest+1;
+for i=2:NTime,
+    TimeGrid(i)=TimeGrid(i-1)+1;
+end
+
+%Bivar = JointPDF(:,1)*PDFTime;
+
+Bivar = zeros(NGrid+1,NTime,NFilter);
+
+for i=1:NGrid+1,
+    for j=1:NTime,
+        for f=1:NFilter,
+            for k=1:NEstimates,
+                %Bivar(i,j) = Bivar(i,j) + exp(-(SCCgrid(i)-SCC(k))/SampleStDev(1)).*exp(-exp(-(SCCgrid(i)-SCC(k))/SampleStDev(1)))/SampleStDev(1) * Weight(k,1) * exp(-0.5*((TimeGrid(j)-Year(k))/StDevYear)^2);
+                if Year(k) >= TimeGrid(j),
+                   Bivar(i,j,f) = Bivar(i,j,f) + exp(-(SCCgrid(i)-SCC(k))/SampleStDev(f)).*exp(-exp(-(SCCgrid(i)-SCC(k))/SampleStDev(f)))/SampleStDev(f) * Weight(k,f) * exp(-(Year(k)-TimeGrid(j)));
+                end
+            end
+        end
+    end
+end
+for f = 1:NFilter,
+    vsum= sum(sum(Bivar(:,:,f)));
+    Bivar(:,:,f) = Bivar(:,:,f)/vsum;
+    MargTime(:,f) = sum(Bivar(:,:,f),1);
+    MargSCC(:,f) = sum(Bivar(:,:,f),2);
+end
+
 clear v*
diff --git a/ConstructCondPDF.m b/ConstructCondPDF.m
@@ -1,27 +1,27 @@
-%% Construct Conditional PDF
-%
-% First version: Richard Tol, 30 October 2012
-% This version: Richard Tol, 30 October 2012
-
-display('Construct conditional PDF');
-
-for i=1:NTime,
-    for f=1:NFilter,
-        if MargTime(i,f) > 0,
-            CondPDF(:,i,f) = Bivar(:,i,f)/MargTime(i,f);
-        else
-            CondPDF(:,i,f) = 0;
-        end
-    end
-end;
-
-CondCDF = CondPDF;
-for i=2:NGrid+1,
-    for j=1:NTime,
-        for f=1:NFilter,
-            CondCDF(i,j,f) = CondCDF(i-1,j,f) + CondPDF(i,j,f);
-        end
-    end
-end
-
+%% Construct Conditional PDF
+%
+% First version: Richard Tol, 30 October 2012
+% This version: Richard Tol, 30 October 2012
+
+display('Construct conditional PDF');
+
+for i=1:NTime,
+    for f=1:NFilter,
+        if MargTime(i,f) > 0,
+            CondPDF(:,i,f) = Bivar(:,i,f)/MargTime(i,f);
+        else
+            CondPDF(:,i,f) = 0;
+        end
+    end
+end;
+
+CondCDF = CondPDF;
+for i=2:NGrid+1,
+    for j=1:NTime,
+        for f=1:NFilter,
+            CondCDF(i,j,f) = CondCDF(i-1,j,f) + CondPDF(i,j,f);
+        end
+    end
+end
+
 clear v*
diff --git a/ConstructHistogram.m b/ConstructHistogram.m
@@ -6,11 +6,14 @@
 
 %%
 hold on
+figure %not used
 plot(SCCgrid(301:2101),JointPDFSilverman(301:2101,1),SCCgrid(301:2101),JointPDFGauss(301:2101,1),SCCgrid(301:2101),JointPDFJSUS(301:2101,1),SCCgrid(301:2101),JointPDFJSU(301:2101,1),SCCgrid(301:2101),JointPDFGaussGauss(301:2101,1),SCCgrid(301:2101),JointPDFGumbelGauss(301:2101,1),SCCgrid(301:2101),JointPDFGumbelWeibull(301:2101,1))
 bar(SCCgrid(301:2101),histo(301:2101,1))
 legend('Normal, Silverman', 'Normal', 'Johnson SU, Silverman', 'Johnson SU', 'Normal, Normal','Gumbel, Normal', 'Gumbel, Weibull', 'Histogram');
 xlabel('dollar per tonne of carbon')
-ylabel('Probability density')
+ylabel('probability density')
+
+hold off
 
 %%
 MISE(1) = histo'*((histo - JointPDFSilverman).*(histo - JointPDFSilverman));
@@ -49,7 +52,7 @@
 mean(7) = SCCgrid*JointPDFGumbelWeibull;
 mean(8) = SCCgrid*histo;
 
-%%
+%% Table S3
 kernelchar = [mean; p0; pl; MISE*1000];
 columnLabels = {'Normal, Silverman', 'Normal', 'Johnson SU, Silverman', 'Johnson SU', 'Normal, Normal','Gumbel, Normal', 'Gumbel, Weibull', 'Observations'};
 rowLabels = {'Average','P(SCC<0)','P(SCC>1186)','MISE'};