% Figures 1 and 2
% August 30, 2015, MWW
% Figures 1 and 2
% ----------- 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 
  h_year = 25;                    % Horizon in years

% ----------- Step 3: Read in data used in analysis      
  % ---- Read in Dataset ----
  vic_data_calendar;
  calvec_save = calvec;
  sfirst = 1947;       % First obs to use (truncated if later)
  slast = 2014.99;     % Last Obs to use (truncated if earlier)
  nper = 4;            % Number of periods per year
  h_per = h_year*nper; % Forecast horizon in periods
   
  % Series names and some parameters 
  sname = 'd_pcpi';
  d_plot = [0.0;0.4;0.7;1.0];
  x_plotmin = -12;
  x_plotmax = 15;
  
%   sname = 'd_gdp_pop';
%   d_plot = [0.0;-0.4;0.2;0.5];
%   x_plotmin = -4;
%   x_plotmax = 6;
    
   % Retrieve data; data series is u, calendar is calvec
      ustr = sname;
      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);  
      h_frac = h_per/dnobs;     % Forecast horizon in fractions of sample
 
      % Output 
      ustr_q = [ustr '_' int2str(q)]; 
      ustr_q = [ustr_q '_matlab'];
      str_out= [outdir ustr_q '_Figures_1_2.out'];
      fid_out = fopen(str_out,'w');   
      
   % ----- Step 4: Print out initial results in output file
    % Plot Series and Projection onto cosine transform
    figure;
        % Figure 1 (i)
        subplot(1,2,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)]);
        plot(calvec,u,'- b','LineWidth',1.0);
        hold on;  
          plot(calvec,uhat,'- k','LineWidth',3.0);
        hold off;
        xlabel('Year')
        
        % Plot Cosine transforms
        %figure 1 (ii)
        subplot(1,2,2);
        newplot;
        tnd = [1:1:q]';
        scatter(tnd,xp1,150,'fill');
        
        set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 22 08]);
        save_g_str = [figdir ustr_q '_fig_1.png'];
        print('-dpng',save_g_str, '-r300');       
        
        
      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);
    
      
      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));
      
        
       % Compute Log-Likelihood for different value of d 
       d_grid = [-0.4:.05:1]';
       n_d = size(d_grid,1);
       gdist = ones(n_d,1)/n_d;  % flat prior
       % Compute Covariance Matrices and save values
       sigma_save = NaN*zeros(q+1,q+1,n_d);
       b = 0;
       c = 0;
       for i = 1:n_d
         d = d_grid(i);
         [sigma] = Sigma_Compute(b,c,d,q,h_frac);
         sigma_save(:,:,i)=sigma;
       end;
       
       
       
       % Carry Out Analyis
       xbar = um;
       X = xp1;
       Xs = X/sqrt(X'*X);
       
       % Construct log-likelihood values
       f_x = NaN*zeros(n_d,1);
       xs_qf = NaN*zeros(n_d,1);
       for i=1:n_d;
        sigma = sigma_save(1:q,1:q,i);
        sigma_inv = inv(sigma);
        f_x(i) = den_invariant(Xs,sigma_inv);
        xs_qf(i) = Xs'*sigma_inv*Xs;
       end;
       f_marg_x = gdist'*f_x;            % Marginal Likelihood
       post_d = (gdist.*f_x)/f_marg_x;   % Posterior
       like_d = f_x;                     % Likelihood values
       llf = log(like_d);                % Log Likihood;
       ii_0 = d_grid == 0;
       llf_0 = llf(ii_0 == 1);
       llf_norm = llf-llf_0;
       
       [tmp,ii] = max(llf_norm);
       fprintf(fid_out,'   MLE of d and value of llf: %6.2f',d_grid(ii));
       fprintf(fid_out,' %6.2f \n',llf_norm(ii));
     
       % Plot Log Likelihood
       figure;
       newplot
       
       subplot(1,2,2);
       % Figure 2: (ii)
       plot(d_grid,llf_norm,'- b','LineWidth',2);
       xlabel(' - d - ')    
       
       % Compute Predictive Density for all values d is dgrid
        ngrid = 1001;                             % Number of grid points for predictive density
        f_yx_save = NaN*zeros(ngrid,n_d);   % Joint density for each horizon
        smin = -40/sqrt(q);   % Min and Max values of Ys; 
        smax = 40/sqrt(q);
        dd = (smax-smin)/(ngrid-1);
        Ys_grid = (smin:dd:smax)';   % grid of values for Ys
        xbarh_grid = xbar + sqrt(X'*X)*Ys_grid;
        Xs_Ys_grid = [repmat(Xs,1,size(Ys_grid,1)); Ys_grid'];   % grid of values for Xs and Ys (Xs held at sample values)
        f_yx = NaN*zeros(size(Ys_grid,1),n_d);  % for joint density
        for i = 1:n_d; 
            sigma = sigma_save(:,:,i); 
            sigma_inv = inv(sigma);
            f_yx(:,i) = den_invariant(Xs_Ys_grid,sigma_inv);
         end;
         % Compute Conditional Density
         f_y_x = f_yx./repmat(f_x',size(f_yx,1),1);  % Conditional density for each value of theta and Ys in grid
         f_y_x_g = f_y_x*post_d;     % Predictive Density using g-prior for year Ys in grid
         
         subplot(1,2,1);
         newplot;
         axis([x_plotmin,x_plotmax,0,inf]);
         box on;
         hold on;
              d = d_plot(1);
              [tmp,ii] = min(abs(d_grid-d));
              f_y_x_d = f_y_x(:,ii);
              plot(xbarh_grid,f_y_x_d,'- k','LineWidth',1);
              d = d_plot(2);
              [tmp,ii] = min(abs(d_grid-d));
              f_y_x_d = f_y_x(:,ii);
              plot(xbarh_grid,f_y_x_d,'-- b','LineWidth',1); 
              d = d_plot(3);
              [tmp,ii] = min(abs(d_grid-d));
              f_y_x_d = f_y_x(:,ii);
              plot(xbarh_grid,f_y_x_d,': k','LineWidth',1);
              d = d_plot(4);
              [tmp,ii] = min(abs(d_grid-d));
              f_y_x_d = f_y_x(:,ii);
              plot(xbarh_grid,f_y_x_d,'-. b','LineWidth',1);
              plot(xbarh_grid,f_y_x_g,'- r','LineWidth',2);
              if strcmp(sname,'d_gdp_pop') == 1;
                legend('d = 0','d = -0.4','d = 0.2','d = 0.5','Bayes','Location','northwest');
              end;
              if strcmp(sname,'d_pcpi') == 1;
                legend('d = 0','d = 0.4','d = 0.7','d = 1.0','Bayes','Location','northwest');
              end;
         hold off;
        
        set(gcf,'PaperUnits','centimeters','PaperPosition',[0 0 22 08]);
        save_g_str = [figdir ustr_q '_fig_2.png'];
        print('-dpng',save_g_str, '-r300');
         
         
      
      type(str_out);
   