general update

CDFtests.m

@@ -0,0 +1,64 @@
+sqn = sqrt(length(SCCgrid));
+for i=18:23,
+    diff = JointCDF(:,1)-JointCDF(:,i);
+    ffid = JointCDF(:,i)-JointCDF(:,1);
+    KS(i-17) = kolmcdf(sqn*max(abs(diff)));
+    K(i-17)  = kolmcdf(sqn*(max(diff)+max(ffid))/2);
+end
+
+%% Filter 18-23 are periods, 1 = null, regardless of discount rate
+res = [5,10,20,50,100];
+for j=1:5,
+    for i=18:23,
+        [ks KC(i-17,j)] = ksdtest(JointCDF(:,i),JointCDF(:,1),res(j));
+    end
+end
+
+columnLabels = {'5', '10', '20','50', '100'};
+rowLabels = {'1982-1995','1996-2001','2002-2006','2007-2013','2014-2017','2018-2021'};
+matrix2latex(KC, 'ks.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');
+
+
+% Filter 24-29 are periods, PRTP = 3, 2 = null
+for j=1:5,
+    for i=24:29,
+        [ks KC3(i-23,j)] = ksdtest(JointCDF(:,i),JointCDF(:,2),res(j));
+    end
+end
+
+columnLabels = {'5', '10', '20','50', '100'};
+rowLabels = {'1982-1995','1996-2001','2002-2006','2007-2013','2014-2017','2018-2021'};
+matrix2latex(KC3, 'ks3.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');
+
+% Filter 30-35 are periods, PRTP = 2, 3 = null
+for j=1:5,
+    for i=30:35,
+        [ks KC2(i-29,j)] = ksdtest(JointCDF(:,i),JointCDF(:,3),res(j));
+    end
+end
+
+columnLabels = {'5', '10', '20','50', '100'};
+rowLabels = {'1982-1995','1996-2001','2002-2006','2007-2013','2014-2017','2018-2021'};
+matrix2latex(KC2, 'ks2.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');
+
+% Filter 36-41 are periods, PRTP = 1, 5 = null
+for j=1:5,
+    for i=36:41,
+        [ks KC1(i-35,j)] = ksdtest(JointCDF(:,i),JointCDF(:,5),res(j));
+    end
+end
+
+columnLabels = {'5', '10', '20','50', '100'};
+rowLabels = {'1982-1995','1996-2001','2002-2006','2007-2013','2014-2017','2018-2021'};
+matrix2latex(KC1, 'ks1.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');
+
+% Filter 42-47 are periods, PRTP = 0, 7 = null
+for j=1:5,
+    for i=42:47,
+        [ks KC0(i-41,j)] = ksdtest(JointCDF(:,i),JointCDF(:,7),res(j));
+    end
+end
+
+columnLabels = {'5', '10', '20','50', '100'};
+rowLabels = {'1982-1995','1996-2001','2002-2006','2007-2013','2014-2017','2018-2021'};
+matrix2latex(KC0, 'ks0.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');

ConstructHistogram.m

@@ -0,0 +1,56 @@
+histo(1) = 0;
+for i = 2:length(SCCgrid)
+    histo(i) = sum(TotalWeight(SCC>=SCCgrid(i-1) & SCC<SCCgrid(i)));
+end
+histo = histo'/sum(histo);
+
+%%
+hold on
+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')
+
+%%
+MISE(1) = histo'*((histo - JointPDFSilverman).*(histo - JointPDFSilverman));
+MISE(2) = histo'*((histo - JointPDFGauss).*(histo - JointPDFGauss));
+MISE(3) = histo'*((histo - JointPDFJSUS).*(histo - JointPDFJSUS));
+MISE(4) = histo'*((histo - JointPDFJSU).*(histo - JointPDFJSU));
+MISE(5) = histo'*((histo - JointPDFGaussGauss).*(histo - JointPDFGaussGauss));
+MISE(6) = histo'*((histo - JointPDFGumbelGauss).*(histo - JointPDFGumbelGauss)); 
+MISE(7) = histo'*((histo - JointPDFGumbelWeibull).*(histo - JointPDFGumbelWeibull));
+MISE(8) = 0;
+
+p0(1) = sum(JointPDFSilverman(SCCgrid<0));
+p0(2) = sum(JointPDFGauss(SCCgrid<0));
+p0(3) = sum(JointPDFJSUS(SCCgrid<0));
+p0(4) = sum(JointPDFJSU(SCCgrid<0));
+p0(5) = sum(JointPDFGaussGauss(SCCgrid<0));
+p0(6) = sum(JointPDFGumbelGauss(SCCgrid<0));
+p0(7) = sum(JointPDFGumbelWeibull(SCCgrid<0));
+p0(8) = sum(histo(SCCgrid<0));
+
+pl(1) = sum(JointPDFSilverman(SCCgrid>1186));
+pl(2) = sum(JointPDFGauss(SCCgrid>1186));
+pl(3) = sum(JointPDFJSUS(SCCgrid>1186));
+pl(4) = sum(JointPDFJSU(SCCgrid>1186));
+pl(5) = sum(JointPDFGaussGauss(SCCgrid>1186));
+pl(6) = sum(JointPDFGumbelGauss(SCCgrid>1186));
+pl(7) = sum(JointPDFGumbelWeibull(SCCgrid>1186));
+pl(8) = sum(histo(SCCgrid>1186));
+
+mean(1) = SCCgrid*JointPDFSilverman;
+mean(2) = SCCgrid*JointPDFGauss;
+mean(3) = SCCgrid*JointPDFJSUS;
+mean(4) = SCCgrid*JointPDFJSU;
+mean(5) = SCCgrid*JointPDFGaussGauss;
+mean(6) = SCCgrid*JointPDFGumbelGauss;
+mean(7) = SCCgrid*JointPDFGumbelWeibull;
+mean(8) = SCCgrid*histo;
+
+%%
+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'};
+matrix2latex(kernelchar, 'kernelchar.tex', 'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.4f', 'size', 'normal');

bootstrap.m

@@ -0,0 +1,41 @@
+%turn off graph and count in ConstructPDF before you run this
+
+for i=1:NEstimates,
+    EstID(i) = i;
+end
+EstID = EstID';
+
+NBootstrap = 1000;
+BootstrapPDF = JointPDF(:,1);
+
+BWeight = zeros(NEstimates,NBootstrap);
+for i=1:NBootstrap,
+    BID = datasample(EstID,NEstimates);
+    for j=1:NEstimates,
+        BWeight(j,i) = sum(BID==j);
+    end
+end
+
+NFilterSave = NFilter;
+NFilter = 1;
+
+JointCDFSave = JointCDF;
+JointPDFSave = JointPDF;
+clear JointCDF JointPDF
+
+for i=1:NBootstrap;
+    display(i)
+    TotalWeight = TotalWeightSave.*BWeight(:,i);
+    ConstructPDF;
+    BootstrapPDF = [BootstrapPDF JointPDF];
+end
+
+%%
+mx = prctile(BootstrapPDF',97.5);
+mn = prctile(BootstrapPDF',2.5);
+%plot(SCCgrid,sum(BootstrapPDF(:,2:end),2)/NBootstrap,SCCgrid,mn,SCCgrid,mx);
+jbfill(SCCgrid(301:2101),mx(301:2101),mn(301:2101),'red','red',0)
+hold on
+plot(SCCgrid(301:2101),sum(BootstrapPDF((301:2101),2:end),2)/NBootstrap);
+xlabel('dollar per tonne of carbon')
+ylabel('Probability density')

ksdtest.m

@@ -0,0 +1,26 @@
+function [ks,p] = ksdtest(testcdf,targetcdf,res)
+%[ks,p] = ksdtest(testcdf,targetcdf,res)
+%ksdtest is a Kolmogorov-Smirnov test on two cumulative distribution
+%functions, testcdf and targetcdf
+%
+%the test is performed on quantiles, defined by res
+%for example, if res=5, the KS test is for the 10th, 30th, 50th, 70th and
+%90th percentiles
+%
+%ks is the test-statistic, p its p-value
+%
+%first version, Richard Tol, 25 April 2021
+%this version, Richard Tol, 8 May 2021
+
+for j=1:res,
+    null(j) = 1/res/2 + (j-1)/res;
+end
+
+for j=1:res,
+    crude(j)=max(testcdf(targetcdf<=null(j)));
+end
+ks = max(abs(crude-null));
+p = 1-kolmcdf(sqrt(res)*ks);
+
+end
+

socialcostcarbon.do

@@ -0,0 +1,40 @@
+regress scc prtp year [iw=cquality]
+regress scc prtp year [iw=cauthor]
+regress scc prtp year [iw=cpaper]
+
+regress scc year [iw=cquality]
+regress scc year [iw=cauthor]
+regress scc year [iw=cpaper]
+
+
+qreg scc prtp year [iw=cquality], q(0.1)
+qreg scc prtp year [iw=cquality], q(0.3)
+qreg scc prtp year [iw=cquality], q(0.5)
+qreg scc prtp year [iw=cquality], q(0.7)
+qreg scc prtp year [iw=cquality], q(0.9)
+qreg scc prtp year [iw=cauthor], q(0.1)
+qreg scc prtp year [iw=cauthor], q(0.3)
+qreg scc prtp year [iw=cauthor], q(0.5)
+qreg scc prtp year [iw=cauthor], q(0.7)
+qreg scc prtp year [iw=cauthor], q(0.9)
+qreg scc prtp year [iw=cpaper], q(0.1)
+qreg scc prtp year [iw=cpaper], q(0.3)
+qreg scc prtp year [iw=cpaper], q(0.5)
+qreg scc prtp year [iw=cpaper], q(0.7)
+qreg scc prtp year [iw=cpaper], q(0.9)
+
+qreg scc year [iw=cquality], q(0.1)
+qreg scc year [iw=cquality], q(0.3)
+qreg scc year [iw=cquality], q(0.5)
+qreg scc year [iw=cquality], q(0.7)
+qreg scc year [iw=cquality], q(0.9)
+qreg scc year [iw=cauthor], q(0.1)
+qreg scc year [iw=cauthor], q(0.3)
+qreg scc year [iw=cauthor], q(0.5)
+qreg scc year [iw=cauthor], q(0.7)
+qreg scc year [iw=cauthor], q(0.9)
+qreg scc year [iw=cpaper], q(0.1)
+qreg scc year [iw=cpaper], q(0.3)
+qreg scc year [iw=cpaper], q(0.5)
+qreg scc year [iw=cpaper], q(0.7)
+qreg scc year [iw=cpaper], q(0.9)