function glogbal_warming_filter
%
% function to use a Kalman Filter to estimate the signal in a sequence of 
% temperature anomalies.  Two data files are available:
%   (1) Global Land-Ocean Temperature Anomalies (deg C) reative to 
%       1951-1980 average
%       d:\temp\Anomalies_NASA_1880_2014.txt
%       and
%   (2) Contiguous 48 U.S. Surface Air Temperature Anomalies (deg C) 
%       relative to 1951-1980 average.
%       d:\temp\Anomalies_US_Surface_Air_Temp_1880_2014.txt
%
% Both data sets are derived from datasets available from the National 
% Aeronautics and Space Administration (NASA) website at:
% http://data.giss.nasa.gov/gistemp/ under Datasets & Images

%==========================================================================
% Function:  global_warming_filter
%
% Author:
%  Rod Deakin, 
%  BONBEACH, VIC, 3196
%  AUSTRALIA
%  email: randm.deakin@gmail.com
%
% Date:
%  Version 1.0  06 April 2015
%  Version 1.1  01 September 2015
%    Moving Average filter added
%
% Remarks:
%  This function uses a Kalman Filter to estimate the signal in global 
%  warming temperature anomalies.
%
% References:
%  [1] Deakin, R.E., 2015, Least Squares and Kalman Filtering,
%        Lecture Notes, September 2015.
%  [2] Deakin, R.E., 2006, The Kalman Filter and Surveying Applications,
%        Lecture Notes, School of Mathematical and Geospatial Sciences,
%        RMIT University, June 2006, 30 pages.
%  [3] Deakin, R.E., 2006, The Kalman Filter: A Look Behind the Scene,
%        Presented at the Regional Surveying Conference, Mildura, 23-25 
%        June 2006, 12 pages.
%
% Arrays:
%  Anomaly   - vector of temperature anomalies
%  B         - Design matrix, dimensions (m,n)
%  b5        - 5-element vector for Moving Average filter
%  corrn     - array of corrections to State vector, dimensions (n,epochs)
%  H         - coefficient matrix of System Driving Noise, dimensions (n,u)
%  I         - Identity matrix, dimensions (n,n)
%  K         - Gain matrix, dimensions (n,m)
%  meas      - vector of measurements, dimensions (epochs,1)
%  Q         - cofactor matrix of measurements, dimensions (n,n)
%  Qmm       - cofactor matrix of Secondary Model, dimensions (n,n)
%  Qww       - cofactor matrix of System Driving Noise, dimensions (n,n)
%  Qxx       - cofactor matrix of State vector, dimensions (n,n)
%  std_xhat  - vector of standard deviations of the filtered State
%  T         - Transition matrix, dimensions (n,n)
%  Temp      - vector of global temperatures
%  U         - cofactor update matrix, dimensions (n,n)
%  xhat      - array containing the State vector at each epoch,
%              dimensions (n,epochs)
%  Year      - vector of years
%  y5        - vector of values for a 5-year moving average filter,
%              dimensions(epochs,1) 
%
% Variables:
%  epochs    - number of years
%  i,j,k     - integer counters
%  m         - number of measurements at particular epoch
%  n         - number of parameters in the State vector 'xhat'
%  std_meas  - standard deviation of measurements
%  u         - number of parameters in the System Driving Noise vector
%==========================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read data from text file %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-------------------------------------------------------------------------
% 1. Call the User Interface (UI) to choose the input data file name
% 2. Concatenate strings to give the path and file name of the input file
% 3. Strip off the extension from the file name to give the rootName
% 4. Add extension ".out" to rootName to give the output filename
% 5. Concatenate strings to give the path and file name of the output file
%-------------------------------------------------------------------------
filepath = strcat('d:\temp\','*.txt');
[infilename,inpathname] = uigetfile(filepath);
infilepath = strcat(inpathname,infilename);
rootName   = strtok(infilename,'.');
outfilename = strcat(rootName,'.out');
outfilepath = strcat(inpathname,outfilename);

%----------------------------------------------------------
% 1. Load the data into an array whose name is the rootName
% 2. set fileTemp = rootName
% 3. Copy columns of data into individual arrays
%----------------------------------------------------------
load(infilepath);
fileTemp = eval(rootName);
Year    = fileTemp(:,1);
Anomaly = fileTemp(:,2);

% determine the number of measurement epochs
epochs = length(Year);

% Set vector of Global Surface Temperatures where the Global Surface 
% Temperature = Anomaly + 14.0 and 14.0 is the mean of the Temperatures for 
% the period 1951-1980.
Temp = Anomaly + 14.0;  % NASA dataset

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the Kalman filter matrices %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 1;  % number of elements in State vector
m = 1;  % number of measurements at each epoch
u = 1;  % number of elements in System Driving Noise vector
% Initialize the State vector.
xhat = zeros(n,epochs);
% Initialize the corrections to State vector.
corrn = zeros(n,epochs);
% Set the State Transition matrix T
T = eye(n);
T(1,1) = 1.0;
% Set the cofactor matrix of the System Driving Noise Qww
Qww = zeros(u,u);
Qww(1,1) = 0.01;
%Qww(1,1) = 0.1;
% Set the coefficient matrix of System Driving Noise H
H = zeros(n,u);
H(1,1) = 1.0;
% Compute cofactor matrix of Dynamic Model by propagation of variances
Qmm = H*Qww*H';
% Initialize vector of measurements
meas = zeros(epochs,1);
% Set the measurement vector to either Temperatures or Temp Anomalies
%meas = Temp;
meas = Anomaly;
% Set the standard deviation of measurements
std_meas = sqrt(0.5);
%std_meas = sqrt(0.2475);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Cofactor matrices and State vector to initial values %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the cofactor matrix of the measurements Q
Q = zeros(m,m);
Q(1,1)  = std_meas^2;
% Set the State cofactor matrix Qxx.  This will be the filtered 
% State cofactor matrix at epoch t2.
Qxx = zeros(n,n);
Qxx(1,1) = std_meas^2;  
% Set the starting estimate of the State vector for Epoch 1
% This will be the filtered Sate at epoch t2
xhat(1,1) = meas(1,1);
% Initialize the Gain matrix
K = zeros(n,m);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start the Kalman Filter at epoch t2 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 2:epochs
  % Compute the predicted State vector.
  xhat(:,k) = T*xhat(:,k-1);
  % Compute the predicted State cofactor matrix
  Qxx = T*Qxx*T' + Qmm;
  % Set the elements of the design matrix B
  B = zeros(m,n);
  B(1,1) = -1;
  % Compute the Gain matrix K
  K = Qxx*B'/(Q + B*Qxx*B');
  % Compute corrections to predicted State
  corrn(:,k) = K*(-meas(k) - B*xhat(:,k));
  % Compute the filtered State vector
  xhat(:,k) = xhat(:,k) + corrn(:,k);
  % Compute the cofactor update matrix
  I = eye(n);
  U = I - K*B;
  % Compute the filtered State cofactor matrix
  Qxx = U*Qxx;
  std_xhat(k) = sqrt(Qxx(1,1));
  
  % Print the State vector, corrections and filtered State cofactor matrix
  fprintf('\n\nEpoch = %3d, Year = %4d, measurement = %8.4f',k,Year(k),meas(k));
  fprintf('\nFiltered State  Corrn         Filtered State cofactor matrix Qxx');
  for i=1:u
      fprintf('\n  %9.4f %10.4f      ',xhat(i,k),corrn(i,k));
      for j=1:u
          fprintf('%12.9f',Qxx(i,j));
      end  
  end  
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% End of Kalman filter %%
%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n\n');

%%%%%%%%%%%%%%
% plot data %%
%%%%%%%%%%%%%%

% Figure 1: Plot of Temperature Anomalies
figure(1);
clf(1);
hold on;
grid on;
box on;
% set scaling for the x- and y-axes
axis([1870 2020 -1.5 1.5]) 
%axis square;
% plot the data
plot(Year,Anomaly,'bo','MarkerSize',3);
% anotate the plot
%title('Global Land-Ocean Temperature Anomalies 1880-2014') 
title('Contiguous 48 U.S. Surface Air Temperature Anomalies 1880-2014') 
xlabel('Year');
ylabel('Temp Anomaly (deg C) relative to 1951-1980 average');

% Figure 2: create a figure window with two plots.  The upper plot shows
% the measurements and the filtered state.  The lower plot shows the
% standard deviation of the filtered state
figure(2);
clf(2);
subplot(2,1,1);
hold on;
grid on;
box on;
plot(meas,'k.','MarkerSize',5);
plot(xhat,'k');
title('Kalman Filter estimate');
xlabel('epoch');
ylabel('state [deg C]');

subplot(2,1,2);
hold on;
grid on;
box on;
plot(std_xhat,'k');
title('Standard deviation of estimate');
xlabel('epoch');
ylabel('std [deg C]');

% Figure 3: Plot of Kalman Filter estimate
figure(3);
clf(3);
hold on;
grid on;
box on;

% set scaling for the x- and y-axes
axis([1870 2020 -1.5 1.5]) 
% plot the data
plot(Year,Anomaly,'bo','MarkerSize',3);
plot(Year,xhat,'k');
% anotate the plot
%title('Global Land-Ocean Temperature Anomalies 1880-2014') 
title({'Contiguous 48 U.S. Surface Air Temperature Anomalies 1880-2014';'with signal (Kalman Filter)'}); 
xlabel('Year');
ylabel('Temp Anomaly (deg C) relative to 1951-1980 average');

% Figure 4: Plot of standard deviation of Kalman Filter estimate
figure(4);
clf(4);
hold on;
grid on;
plot(std_xhat,'k');
hold off;
title('Standard deviation of estimate');
xlabel('epoch');
ylabel('std [deg C]');

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Moving Average Filter %%
%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5-year moving average
b5  = [1/5 1/5 1/5 1/5 1/5];
y5  = filter(b5,1,Anomaly);

% centre the averages about mid-values
for k = 3:epochs-2
    y5(k) = y5(k+2);     
end
y5(1) = 0;
y5(2) = 0;
y5(epochs-1) = 0;
y5(epochs)   = 0;

% Figure 5: Moving average filter
figure(5);
clf(5);
hold on;
grid on;
box on;
% set scaling for the x- and y-axes
axis([1870 2020 -1.5 1.5]) 
% plot the data
plot(Year,Anomaly,'bo','MarkerSize',3);
% plot centred 5-year moving average
plot(Year(3:epochs-2),y5(3:epochs-2),'k-','LineWidth',1);
% anotate the plot
title({'Contiguous 48 U.S. Surface Air Temperature Anomalies 1880-2014';'with centred 5-Year Moving Average'}); 
xlabel('Year');
ylabel('Temp Anomaly (deg C) relative to 1951-1980 average');


