% lr_main_quarterly_returns_1947.m
% August 31, 2015, MWW
% This program carries on long-run prediction Results for Quarterly Data
% ----------- Step 1: Set up parameters and so forth for run
clear all;
small = 1.0e-10;
  
  % -- File Directories 
  outdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/out/';
  figdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/fig/';
  matdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/mat/';;
  mat_save_dir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/mat_save/';  % Location for Saving Various Files
 
  
  % -- Parameters
  q = 12;                         % value of q to use 
  clevel_vec = [0.67;0.90];  % Confidence Level
  hyears_tab = [10;25;50;75];     % Horizons to Tabulate
  
  % -- Parameters to govern calculation of covariance matrices, grid of BCD values and so forth
  i_first = 0;  % = 1, the first time this program is used.  In this case, the program sets up various matrices
                %      (cov. matrices, BCD values, horizon fractions, etc.) that can be used for future calls.

                
  % ---------- Step 2: Set up BCD parameters, Sigma matrices and so forth
  if i_first == 1;
    lr_main_setup_matrices;
  end;                 

% ----------- Step 3: Read in data used in analysis      
  % ---- Read in Dataset ----
  vic_data_calendar_returns;
  calvec_save = calvec;
  sfirst = 1926;  % First obs to use (truncated if later)
  slast = 2014.99;  % Last Obs to use (truncated if earlier)
  nper = 4;         % Number of periods per year

  % list of series to use
    slist = {...
             'vwret_q_real' ... 
             }'; 
 
         
 % --------------------- Begin Loop over Series -----------
  for iseries = 1:size(slist,1);
  
   % Retrieve data; data series is u, calendar is calvec
      ustr = char(slist(iseries));
      eval(['u = ' ustr ';']);
      tmp = packr([calvec_save,u]);
      calvec = tmp(:,1);
      u = tmp(:,2);
   % Find Series over sample period 
      ismpl = (calvec >= sfirst)&(calvec <= slast);
      calvec = calvec(ismpl==1,:);
      u = u(ismpl==1,:);
      dnobs = size(calvec,1);  
 
   % ----- Step 4: Print out initial results in output file
     % Some Summary Statistics 
      [up,xp1] = xp_compute(u,q);
      nma = floor(1.5*(dnobs^.333)+1);
      [um,var_um]=hacm(u,nma,1);
      lrcov = var_um*dnobs;          % Newey-West Var
      var_lr = dnobs*mean(xp1.^2);
    % Output 
      ustr_q = [ustr '_' int2str(q) '_' int2str(sfirst) '_' int2str(slast) '_']; 
      str_out= [outdir ustr_q '_lr_main.out'];
      fid_out = fopen(str_out,'w');
      fprintf(fid_out,[' Date and Time: ' datestr(now) '\n']);
      fprintf(fid_out,[' Series name: ' ustr '\n']);
      fprintf(fid_out,'   First obs = %6.2f\n',calvec(1));
      fprintf(fid_out,'   Last obs = %6.2f\n',calvec(size(calvec,1)));
      fprintf(fid_out,'   Number of Obs %4.0f\n',dnobs);
      fprintf(fid_out,'   Number of Years %4.1f\n',dnobs/nper);
      fprintf(fid_out,'   Number of Cosine Transformations %4d\n',q);
      fprintf(fid_out,'   Cutoff periodicity (in years) = %4.1f\n',(2*dnobs/q)/nper);
      fprintf(fid_out,'   Sample Mean = %6.2f\n',um);
      fprintf(fid_out,'   Sample Std. Dev. = %6.2f\n',std(u));
      fprintf(fid_out,'   LR SD (New-West) = %6.2f\n',sqrt(lrcov));
      fprintf(fid_out,'   LR SD (Low-Freq I(0)) = %6.2f\n',sqrt(var_lr));
   % ---- Step 5: Compute I(0) Prediction Sets
      fprintf(fid_out,'   I(0) Prediction Sets for 25-year horizon \n');
      r = (nper*25)/dnobs;  % horizon fraction
      for ic = 1:size(clevel_vec,1);
          tq = tinv((1-(1-clevel_vec(ic))/2),q);
          del = tq*sqrt((1+(1/r))*(1/dnobs)*var_lr);
          fprintf(fid_out,'      Level %6.2f: (%4.1f,%4.1f)\n',clevel_vec(ic),um-del,um+del);
      end;
   % ---- Step 6: Compute bcd prediction sets 
      %  ----- Read in variables
      str_tmp = 'hfrac_vec'; load([matdir str_tmp '_' int2str(q)]);    
      str_tmp = 'bcd_mat'; load([mat_save_dir str_tmp '_' int2str(q)]);    
      str_tmp = 'gdist'; load([mat_save_dir str_tmp '_' int2str(q)]);  
      str_tmp = 'cvalue_mat'; load([mat_save_dir str_tmp '_' int2str(q)]); 
      str_tmp = 'lfd_mat'; load([mat_save_dir str_tmp '_' int2str(q)]); 
      str_tmp = 'w_mat'; load([mat_save_dir str_tmp '_' int2str(q)]);
      str_tmp = 'sig_matlab_save'; load([mat_save_dir str_tmp '_' int2str(q)]); 
      % --- Some Dimensions
      n_bcd = size(bcd_mat,1);
      n_clevels = size(clevel_vec,1);
      n_horizons = size(hfrac_vec,1);
      
      % Choose sigma matrix to use;
      sigma_save = sig_matlab_save;
      
      % Carry out calculations 
        [Bayes_pred_sets,MN_pred_sets,like_bcd,post_bcd,sigma_save_output] ... 
         = lr_pred_set(u,q,hfrac_vec,clevel_vec,bcd_mat,sigma_save,gdist,w_mat,lfd_mat,cvalue_mat);
         % Save sigma matrix if i_um = 0 and it hasn't been computed before;
         if (size(sig_matlab_save,1)==1);
           sig_matlab_save = sigma_save_output;
           str_tmp = [mat_save_dir 'sig_matlab_save_' int2str(q)];
           save(str_tmp,'sig_matlab_save');
         end;
     
     % --------- Step 7 --- Print out log-likelihood for fractional model at particular values of d
       ii = (bcd_mat(:,1)==0).*(bcd_mat(:,2)==0);
       dvec = bcd_mat(ii==1,3);
       like_dvec = like_bcd(ii==1);
       id0 = (dvec==0);
       rel_ll = log(like_dvec/like_dvec(id0==1));
       fprintf(fid_out,'\n\n   Log Likelihood values for d values \n');
       for id =1:size(dvec,1);
           fprintf(fid_out,'       %3.1f , %4.1f \n',dvec(id),rel_ll(id));
       end;
         
    % ---- Step 8  Print Predictions Sets for different horizons
    
    % -- By horizon fraction (not reported in paper)
    for ih = 1:n_horizons;
      fprintf(fid_out,'%4.2f~',hfrac_vec(ih));
      for ic = 1:n_clevels;
         fprintf(fid_out,'(%4.2f,%4.2f)',Bayes_pred_sets(1:2,ic,ih));
         fprintf(fid_out,'[%4.2f,%4.2f]',MN_pred_sets(1:2,ic,ih));
         if ic == n_clevels;
           fprintf(fid_out,'\n');
         else;
           fprintf(fid_out,'~');
         end;
      end;
    end; 

    
    % ----- Results for Tables 4, 5, 6 -- Calculate Results via linear interpolation
      n_rpt = size(hyears_tab,1);
      Bayes_pred_rpt = NaN*zeros(2,n_clevels,n_rpt);
      MN_pred_rpt = NaN*zeros(2,n_clevels,n_rpt);
      for ic = 1:n_clevels;
        for irpt = 1:n_rpt;
           hy = hyears_tab(irpt);
           hf = (nper*hy)/dnobs;
    	     % Linear interpolation to find results 
 	         % Find largest hfract <= hf 
 	         tmp = (hfrac_vec).*(hfrac_vec <= hf);
 	         if max(tmp) == 0;
  	           ih1 = 1; ih2 = ih1; wght1=1;
             elseif max(tmp) == max(hfrac_vec);
  	           ih1 = size(hfrac_vec,1); ih2 = ih1; wght1 = 1;
             else;
               [junk,ih1] = max(tmp); ih2 = ih1+1; 
                wght1 = 1-((hf-hfrac_vec(ih1))/(hfrac_vec(ih2)-hfrac_vec(ih1)));
             end;
           hf1 = hfrac_vec(ih1)*wght1 + hfrac_vec(ih2)*(1-wght1);
           Bayes_pred_rpt(:,ic,irpt) = wght1*Bayes_pred_sets(:,ic,ih1)+(1-wght1)*Bayes_pred_sets(:,ic,ih2);
           MN_pred_rpt(:,ic,irpt) = wght1*MN_pred_sets(:,ic,ih1)+(1-wght1)*MN_pred_sets(:,ic,ih2);     
        end;
      end;
      % Print Out Results 
      fprintf(fid_out,'\n\n Predictions Sets (Bayes) [MN] \n');
      for ic = 1:n_clevels;
          fprintf(fid_out,'  Coverage level %4.2f \n',clevel_vec(ic));
          fprintf(fid_out,'    %3.0f,',hyears_tab(1:end-1));
          fprintf(fid_out,'    %3.0f \n',hyears_tab(end));
          for irpt = 1:n_rpt;
             fprintf(fid_out,'    (%4.1f,%4.1f)',Bayes_pred_rpt(1:2,ic,irpt));
             if irpt < n_rpt;
                fprintf(fid_out,',');  
              else;
                fprintf(fid_out,'\n');
              end;
          end;
          for irpt = 1:n_rpt;
             fprintf(fid_out,'    [%4.1f,%4.1f]',MN_pred_rpt(1:2,ic,irpt));
             if irpt < n_rpt;
                fprintf(fid_out,',');  
             else;
                fprintf(fid_out,'\n');
             end;
          end;
      end; 

      % ------------- Step 9 Print graph of series and low frequency component for Appendix A.1
        psi = psi_compute(dnobs,q);
        um = mean(u);
        udm = u - um;
        xp1 = psi'*udm; 
        uhat = psi*(inv(psi'*psi))*xp1;
        uhat = uhat + um;
        plot(calvec,u);
        v = axis;
        axis([1940 2020 v(3) v(4)]);
        hold on;
          plot(calvec,u,'- b','LineWidth',1);
          plot(calvec,uhat,'- r','LineWidth',2);
        hold off;
        save_g_str = [figdir ustr_q '_series_uproj.jpg'];
        saveas(gcf,save_g_str);
   
        
        % ------------- Step 9 Print graph of series and low frequency component for Appendix A.1
        psi = psi_compute(dnobs,q);
        um = mean(u);
        udm = u - um;
        xp1 = psi'*udm; 
        uhat = psi*(inv(psi'*psi))*xp1;
        uhat = uhat + um;
      
      % Set up Plot 
      figure;
      
      % Plot Standardized Cosine transforms
          tnd = [1:1:q]';
          Xs = xp1/sqrt(xp1'*xp1);
          box on;
          subplot(2,3,1);
          plot(tnd,Xs,'ok','MarkerSize',10,'MarkerFaceColor','b');  
          title(['(a) Cosine Transform']);
          
          % Plot Log Likelihood
          box on;
          subplot(2,3,2);
          plot(dvec,rel_ll,'-ok','MarkerSize',10,'LineWidth',3);
          title(['(b) Log-Likelihood for d']);
          
      
      % ------------- Step 9 Print graph of series and forecasts for figures in Appendix A.2
          % Step 9.1 ... Add 75 years to calendar @
          calvec2 = (calvec(end)+1/nper:1/nper:calvec(end)+75)';
          calvec_a = [calvec;calvec2];
          % Series to be plotted
          dnobs_a = size(calvec_a,1);
          u_actual = NaN*zeros(dnobs_a,1);
          u_actual(1:dnobs) = u;
          u_lf = NaN*zeros(dnobs_a,1);
          u_lf(1:dnobs) = uhat;
          bayes_ci = NaN*zeros(dnobs_a,2,n_clevels);
          MN_ci = NaN*zeros(dnobs_a,2,n_clevels);
          I0_ci = NaN*zeros(dnobs_a,2,n_clevels);
          for ic = 1:n_clevels;
            for t = ceil(1.1*dnobs):dnobs+75*nper;
              hfrac_act = (t-dnobs)/dnobs;  % rho
 	            % Linear interpolation to find results
              % Find largest hfract <= hf 
 	            tmp = (hfrac_vec).*(hfrac_vec <= hfrac_act);
 	            if max(tmp) == 0;
  	               ih1 = 1; ih2 = ih1; wght1=1;
              elseif max(tmp) == max(hfrac_vec);
  	               ih1 = size(hfrac_vec,1); ih2 = ih1; wght1 = 1;
              else;
                   [junk,ih1] = max(tmp); ih2 = ih1+1; 
                   wght1 = 1-((hfrac_act-hfrac_vec(ih1))/(hfrac_vec(ih2)-hfrac_vec(ih1)));
              end;
              bayes_ci(t,:,ic)=wght1*Bayes_pred_sets(:,ic,ih1)+(1-wght1)*Bayes_pred_sets(:,ic,ih2);
              MN_ci(t,:,ic)=wght1*MN_pred_sets(:,ic,ih1)+(1-wght1)*MN_pred_sets(:,ic,ih2);  
              
              tq = tinv((1-(1-clevel_vec(ic))/2),q);
              del = tq*sqrt((1+(1/hfrac_act))*(1/dnobs)*var_lr);
              I0_ci(t,1,ic) = um-del;
              I0_ci(t,2,ic) = um+del;
          
            end;
          end;       
          
          % Step 9.2 construct plots
          subplot(2,3,4);
          plot(calvec_a,[u_actual MN_ci(:,1,2) MN_ci(:,2,2)]);
          v = axis;
          axis([1945 2015 v(3) v(4)]);
          hold on;
            plot(calvec_a,u_actual,'- b','LineWidth',1);
            plot(calvec_a,u_lf,'- k','LineWidth',2);
          hold off;
          title(['(c) Data and Low-frequency projection']);
          %plot(calvec_a,u_actual);
          
        
          for ic = 1:n_clevels;
              subplot(2,3,4+ic);
              axis([2015 2095 v(3) v(4)]); 
              hold on; 
                 plot(calvec_a,bayes_ci(:,1,ic),'- b','LineWidth',3);
                 plot(calvec_a,bayes_ci(:,2,ic),'- b','LineWidth',3);
                 plot(calvec_a,MN_ci(:,1,ic),'- b','LineWidth',1);
                 plot(calvec_a,MN_ci(:,2,ic),'- b','LineWidth',1);
                 plot(calvec_a,I0_ci(:,1,ic),'-- r','LineWidth',1);
                 plot(calvec_a,I0_ci(:,2,ic),'-- r','LineWidth',1);
              hold off;
              if ic == 1;
                 title(['(d) 67% Prediction Sets']);
              else;
                 title(['(d) 90% Prediction Sets']);
              end;  
              tmp = floor(100*clevel_vec(ic));  
          end;
          
          set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 28 14]);
          save_g_str = [figdir ustr_q '_appendix.png'];
          print('-dpng',save_g_str, '-r300');
                   
      type(str_out);
  end;
       
 