phi = 2.85;
a = 0.72;
d = 2.58;
e = 1e-10;

X = fatdata(:,1:13);
Y = fatdata(:,14);
[N,p] = size(X);

display('============= full data ===================')
          x = X;
          y = Y;
          n = N;
          xx = [ones(n,1),x];
          %handle alasso
          betalasso = alasso(center(x),center(y));
          betalasso = [mean(y)-mean(x)*betalasso;betalasso];
          sig2lasso = (y-xx*betalasso)'*(y-xx*betalasso)/n;
          
 		  % handle plasso	
          obj = adaplasso_est(x,y);
          betaplasso = obj.beta;
          sig2plasso = obj.sig2;
          
          % handle wplasso
          obj = adapwplasso_est(x,y);
          betawplasso = obj.beta;
          sig2wplasso = obj.sig2;
          display('===== estimated coefficients by aLasso, pLasso, wpLasso =====')
          [betalasso,betaplasso,betawplasso]
display('==========================================')            
display('============= case I ===================')
display('==========================================')  
          n = 125;
          ntest = N-n;    
          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 alasso
          
          betalasso = alasso(center(x),center(y));
          betalasso = [mean(y)-mean(x)*betalasso;betalasso];
          sig2lasso = (y-xx*betalasso)'*(y-xx*betalasso)/n;
          PPSlasso = (1/2)*log(sig2lasso)+(xxtest*betalasso-ytest)'*(xxtest*betalasso-ytest)/(2*sig2lasso*ntest);
          
          sighat = sqrt(sig2lasso);
          muv = xxtest*betalasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSalasso = aa/ntest;
          
 		  % handle plasso	
          obj = adaplasso_est(x,y);
          betaplasso = obj.beta;
          sig2plasso = obj.sig2;
          PPSplasso = (1/2)*log(sig2plasso)+(xxtest*betaplasso-ytest)'*(xxtest*betaplasso-ytest)/(2*sig2plasso*ntest);

          sighat = sqrt(sig2plasso);
          muv = xxtest*betaplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSplasso = aa/ntest;
          
          
          % handle wplasso
          obj = adapwplasso_est(x,y);
          betawplasso = obj.beta;
          sig2wplasso = obj.sig2;
          phi = 2.85; a = 0.72; d = 2.58;
          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 = PPS/ntest;
          
          sighat = sqrt(sig2wplasso);
          muv = xxtest*betawplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSwplasso = aa/ntest;

          display('===== estimated coefficients by aLasso, pLasso, wpLasso =====')
          [betalasso,betaplasso,betawplasso]
          display('===== PPS by aLasso, pLasso, wpLasso =====')
          [PPSlasso,PPSplasso,PPSwplasso]
          display('===== CRPS by aLasso, pLasso, wpLasso =====')
          [CRPSalasso,CRPSplasso,CRPSwplasso]
display('==========================================')    
display('============= case II ===================')
display('==========================================')  
          n = 125;
          ntest = N-n;    
          xtest = X(1:ntest,:);
          ytest = Y(1:ntest);
          x = X((ntest+1):N,:);
          y = Y((ntest+1):N);
          xx = [ones(n,1),x];
          xxtest = [ones(ntest,1),xtest];

          %handle alasso
          
          betalasso = alasso(center(x),center(y));
          betalasso = [mean(y)-mean(x)*betalasso;betalasso];
          sig2lasso = (y-xx*betalasso)'*(y-xx*betalasso)/n;
          PPSlasso = (1/2)*log(sig2lasso)+(xxtest*betalasso-ytest)'*(xxtest*betalasso-ytest)/(2*sig2lasso*ntest);
          
          sighat = sqrt(sig2lasso);
          muv = xxtest*betalasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSalasso = aa/ntest;
          
 		  % handle plasso	
          obj = adaplasso_est(x,y);
          betaplasso = obj.beta;
          sig2plasso = obj.sig2;
          PPSplasso = (1/2)*log(sig2plasso)+(xxtest*betaplasso-ytest)'*(xxtest*betaplasso-ytest)/(2*sig2plasso*ntest);

          sighat = sqrt(sig2plasso);
          muv = xxtest*betaplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSplasso = aa/ntest;
          
          
          % handle wplasso
          obj = adapwplasso_est(x,y);
          betawplasso = obj.beta;
          sig2wplasso = obj.sig2;
          phi = 2.85; a = 0.72; d = 2.58;
          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 = PPS/ntest;
          
          sighat = sqrt(sig2wplasso);
          muv = xxtest*betawplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSwplasso = aa/ntest;
          display('===== estimated coefficients by aLasso, pLasso, wpLasso =====')
          [betalasso,betaplasso,betawplasso]
          display('===== PPS by aLasso, pLasso, wpLasso =====')
          [PPSlasso,PPSplasso,PPSwplasso]
          display('===== CRPS by aLasso, pLasso, wpLasso =====')
          [CRPSalasso,CRPSplasso,CRPSwplasso]
display('==========================================')  
display('============= case III ===================')
display('==========================================')  
          v = randperm(N);
          X = fatdata(v,1:13);
          Y = fatdata(v,14);
          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 alasso
          
          betalasso = alasso(center(x),center(y));
          betalasso = [mean(y)-mean(x)*betalasso;betalasso];
          sig2lasso = (y-xx*betalasso)'*(y-xx*betalasso)/n;
          PPSlasso = (1/2)*log(sig2lasso)+(xxtest*betalasso-ytest)'*(xxtest*betalasso-ytest)/(2*sig2lasso*ntest);
          
          sighat = sqrt(sig2lasso);
          muv = xxtest*betalasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSalasso = aa/ntest;
          
 		  % handle plasso	
          obj = adaplasso_est(x,y);
          betaplasso = obj.beta;
          sig2plasso = obj.sig2;
          PPSplasso = (1/2)*log(sig2plasso)+(xxtest*betaplasso-ytest)'*(xxtest*betaplasso-ytest)/(2*sig2plasso*ntest);

          sighat = sqrt(sig2plasso);
          muv = xxtest*betaplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSplasso = aa/ntest;
          
          
          % handle wplasso
          obj = adapwplasso_est(x,y);
          betawplasso = obj.beta;
          sig2wplasso = obj.sig2;
          phi = 2.85; a = 0.72; d = 2.58;
          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 = PPS/ntest;
          
          sighat = sqrt(sig2wplasso);
          muv = xxtest*betawplasso;
          aa = 0;
          for kk = 1:ntest
              aa = aa+crps(muv(kk),sighat,ytest(kk));
          end
          CRPSwplasso = aa/ntest;
          display('===== estimated coefficients by aLasso, pLasso, wpLasso =====')
          [betalasso,betaplasso,betawplasso]
          display('===== PPS by aLasso, pLasso, wpLasso =====')
          [PPSlasso,PPSplasso,PPSwplasso]
          display('===== CRPS by aLasso, pLasso, wpLasso =====')
          [CRPSalasso,CRPSplasso,CRPSwplasso]
