% This script is to compare the (adaptive and non-adaptive) plasso, wplasso 
% to lasso and alasso in terms of PPS, MSE, NZE and CRPS in linear
% regression.
%
% small-p case and normal predictors.
%
% This was used to produce Table 1 in the paper "The predictive Lasso" by
% Tran, Nott and Leng.
%
clear all
warning off 
n = 200; sig = 3;
rept = 20;                 % number of replications

p = 8;
tbeta = [3,1.5,0,0,0.5,.5,0,0]';
Sig = ones(p);
for ii = 1:p
	for jj = 1:p
		if (ii~=jj) Sig(ii,jj) = 0.5^(abs(ii-jj));
        end
    end
end

truebeta = [2;tbeta];
trueflag = abs(truebeta)>0;
phi = 2.85;
a = 0.72;
d = 2.58;
e = 1e-10;

      ntest = n;
      N = n+ntest;
      for irep = 1:rept
          X = mvnrnd(zeros(1,p),Sig,N);
          %X = mvtrnd(Sig,1.5,N);
          Y = 2+X*tbeta + sig*randn(N,1);
          x = X(1:n,:);
          y = Y(1:n);
          xtest = X((n+1):N,:);
          ytest = Y((n+1):N);
          xx = [ones(n,1),x];
          xxtest = [ones(ntest,1),xtest];

          % handle lasso
 
          betahat = lasso(center(x),center(y));
          betahat = [mean(y)-mean(x)*betahat;betahat];
          sig2hat = (y-xx*betahat)'*(y-xx*betahat)/n;

          PPSlasso(irep) = (1/2)*log(sig2hat)+(xxtest*betahat-ytest)'*(xxtest*betahat-ytest)/(2*sig2hat*ntest);
          
          NZElasso(irep) = sum(abs(betahat)<e);
          
          MSElasso(irep) = (truebeta-betahat)'*(truebeta-betahat); 

          sighat = sqrt(sig2hat);
          muv = xxtest*betahat;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSlasso(irep) = aa/ntest;
         
          CFRlasso(irep) = sum((abs(betahat)>e) == trueflag)==p+1;
          
          %handle alasso
          betahat = alasso(center(x),center(y));
          betahat = [mean(y)-mean(x)*betahat;betahat];
          sig2hat = (y-xx*betahat)'*(y-xx*betahat)/n;

          PPSalasso(irep) = (1/2)*log(sig2hat)+(xxtest*betahat-ytest)'*(xxtest*betahat-ytest)/(2*sig2hat*ntest);
          
          NZEalasso(irep) = sum(abs(betahat)<e);
          
          MSEalasso(irep) = (truebeta-betahat)'*(truebeta-betahat); 

          sighat = sqrt(sig2hat);
          muv = xxtest*betahat;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSalasso(irep) = aa/ntest;

          CFRalasso(irep) = sum((abs(betahat)>e) == trueflag)==p+1;
          
 		  % handle non-adaptive plasso	
          obj = naplasso_est(x,y);
          betahat = obj.beta;
          sig2hat = obj.sig2;
          
          PPSnaplasso(irep) = (1/2)*log(sig2hat)+(xxtest*betahat-ytest)'*(xxtest*betahat-ytest)/(2*sig2hat*ntest);
          
          NZEnaplasso(irep) = sum(abs(betahat)<e);
          
          MSEnaplasso(irep) = (truebeta-betahat)'*(truebeta-betahat);
          
          sighat = sqrt(sig2hat);
          muv = xxtest*betahat;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSnaplasso(irep) = aa/ntest;
          
          CFRnaplasso(irep) = sum((abs(betahat)>e) == trueflag)==p+1;
 		  % handle adaptive plasso	
          obj = adaplasso_est(x,y);
          betaplasso = obj.beta;
          sig2plasso = obj.sig2;
          
          PPSaplasso(irep) = (1/2)*log(sig2plasso)+(xxtest*betaplasso-ytest)'*(xxtest*betaplasso-ytest)/(2*sig2plasso*ntest);
          
          NZEaplasso(irep) = sum(abs(betaplasso)<e);
          
          MSEaplasso(irep) = (truebeta-betaplasso)'*(truebeta-betaplasso);
          
          sigplasso = sqrt(sig2plasso);
          muv = xxtest*betaplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sigplasso,ytest(kk));
          end
          CRPSaplasso(irep) = aa/ntest;
          CFRaplasso(irep) = sum((abs(betaplasso)>e) == trueflag)==p+1;
          
          % handle wplasso
          obj = adapwplasso_est(x,y);
          betawplasso = obj.beta;
          sig2wplasso = obj.sig2;

          m = xx'*xx\xx'*y;
          m(2:(p+1)) = 0;
          Vinv = diag([1/var(y),var(x)./phi^2]);
          Vhat = inv(Vinv+xx'*xx);
          wtest = zeros(ntest,1);
          for i = 1:ntest
              wtest(i) = 1+xxtest(i,:)*Vhat*xxtest(i,:)';
          end
          PPS = 0;
          for i = 1:ntest
              PPS = PPS + (1/2)*log(sig2wplasso*wtest(i))+(1/2)*(ytest(i)-xxtest(i,:)*betawplasso)^2/(sig2wplasso*wtest(i));
          end
          PPSwplasso(irep) = PPS/ntest;

          NZEwplasso(irep) = sum(abs(betawplasso)<e);
 
          MSEwplasso(irep) = (truebeta-betawplasso)'*(truebeta-betawplasso);
          
          sigwplasso = sqrt(sig2wplasso);
          muv = xxtest*betawplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sigwplasso,ytest(kk));
          end
          CRPSwplasso(irep) = aa/ntest;
          CFRwplasso(irep) = sum((abs(betawplasso)>e) == trueflag)==p+1;
      end
      
      fprintf('============= Lasso - aLasso - non-adaptive pLasso - adaptive pLasso - wpLasso===============\n');
      fprintf('    *n = %1d, sigma = %1.0f*\n',n, sig);
      fprintf('******** PPS: mean and sd **********\n');
      [mean(PPSlasso),sqrt(var(PPSlasso)),mean(PPSalasso),sqrt(var(PPSalasso)),...
          mean(PPSnaplasso),sqrt(var(PPSnaplasso)),mean(PPSaplasso),sqrt(var(PPSaplasso)),mean(PPSwplasso),sqrt(var(PPSwplasso))]
      fprintf('******** MSE: mean and sd **********\n');
      [mean(MSElasso),sqrt(var(MSElasso)),mean(MSEalasso),sqrt(var(MSEalasso)),...
          mean(MSEnaplasso),sqrt(var(MSEnaplasso)),mean(MSEaplasso),sqrt(var(MSEaplasso)),mean(MSEwplasso),sqrt(var(MSEwplasso))]
      fprintf('******** NZE: mean and sd **********\n');
      [mean(NZElasso),sqrt(var(NZElasso)),mean(NZEalasso),sqrt(var(NZEalasso)),...
          mean(NZEnaplasso),sqrt(var(NZEnaplasso)),mean(NZEaplasso),sqrt(var(NZEaplasso)),mean(NZEwplasso),sqrt(var(NZEwplasso))]
      fprintf('******** CPRS: mean and sd *********\n');
      [mean(CRPSlasso),sqrt(var(CRPSlasso)),mean(CRPSalasso),sqrt(var(CRPSalasso)),...
          mean(CRPSnaplasso),sqrt(var(CRPSnaplasso)),mean(CRPSaplasso),sqrt(var(CRPSaplasso)),mean(CRPSwplasso),sqrt(var(CRPSwplasso))]
      fprintf('******** CFR *********\n');
      [mean(CFRlasso),mean(CFRalasso),mean(CFRnaplasso),mean(CFRaplasso),mean(CFRwplasso)]

% boxplots
PPS = [PPSlasso',PPSalasso',PPSnaplasso',PPSaplasso',PPSaplasso'];
MSE = [MSElasso',MSEalasso',MSEnaplasso',MSEaplasso',MSEaplasso'];
CRPS = [CRPSlasso',CRPSalasso',CRPSnaplasso',CRPSaplasso',CRPSaplasso'];
NZE = [NZElasso',NZEalasso',NZEnaplasso',NZEaplasso',NZEaplasso'];
subplot(2,2,1);
boxplot(PPS)
title('PPS')
xlabel('')
set(gca,'XTickLabel',{'L','aL','npL','pL','wpL'})

subplot(2,2,2);
boxplot(MSE)
title('MSE')
xlabel('')
set(gca,'XTickLabel',{'L','aL','npL','pL','wpL'})

subplot(2,2,3);
boxplot(CRPS)
title('CRPS')
xlabel('')
set(gca,'XTickLabel',{'L','aL','npL','pL','wpL'})

subplot(2,2,4);
boxplot(NZE)
title('NZE')
xlabel('')
set(gca,'XTickLabel',{'L','aL','npL','pL','wpL'})
