function a=warmingworld
% a=WARMINGWORLD
%
% Loads and plots data ultimately from NASA, proximally from
% https://www.bloomberg.com/graphics/2015-whats-warming-the-world/
%
% Needs updating to take input on which variables to consider
%
% Last modified by fjsimons-at-alum.mit.edu, 01/11/2023

% Define the directory where you keep the data
pname='/Users/fjsimons/Downloads';

% Load the data tables
o=readtable(fullfile(pname,'observed.csv'));
f=readtable(fullfile(pname,'forcings.csv'));

% Work from 1880 onwards to 2005
loy=1880; hiy=2005;
selekto=o.Year>=loy & o.Year<=hiy;
selektf=f.Year>=loy & f.Year<=hiy;
oyear=o.Year(selekto);
fyear=f.Year(selektf);

% Pick out columns and standardize
oamean=[o.Annual_Mean(selekto)-mean(o.Annual_Mean(selekto))]/std(o.Annual_Mean(selekto));
fghg=[f.GreenhouseGases(selektf)-mean(f.GreenhouseGases(selektf))]/std(f.GreenhouseGases(selektf));
fsolar=[f.Solar(selektf)-mean(f.Solar(selektf))]/std(f.Solar(selektf));
fvolc=[f.Volcanic(selektf)-mean(f.Volcanic(selektf))]/std(f.Volcanic(selektf));

% Plot a figure
clf
plot(oyear,oamean,'k','LineWidth',2)
hold on
plot(fyear,fghg)
plot(fyear,fsolar)
plot(fyear,fvolc)
hold off

% Add a legend
l=legend('Temperature','Greenhouse','Solar','Volcanic',...
         'Location','NorthWest');

% Now perform a multiple linear regression for the model
% oamean=a*fghg+b*fsolar+c*fvolc;
% We call the unknown parameter vector [a b c] here pvars
G=[fsolar fvolc];

% Perform the inversion for the unkown multipliers
pvars=pinv(G)*oamean;

% Make a PREDICTION based on these multipliers
preds=G*pvars;

% And plot that right on top of the observations
hold on
plot(oyear,preds,'b','LineW',1.5)
hold off

% Label properly
grid on
xlim([loy hiy])
xlabel('year')
ylabel('Temperature | Greenhouse | Solar | Volcanic')

% Calculate the rms residuals
res=oamean-preds;

% Make histograms of the residuals
% You do this

% Plot the correlation coefficient of the predictions with the
% observations
[R,P]=corrcoef(preds,oamean);
% Pull out the correlation coefficients
R=R(2);
% Pull out the probabilities of these correlations, or higher, arising by
% chance when in fact, there are no correlations
P=P(2);

% Data variance (Total Variance)
dvar=var(oamean);
% Variance of the residuals (Unexplained Variance)
rvar=var(res);
% Variance reduction is the proportion of Explained Variance compared to
% the Total Variance
vardux=(dvar-rvar)/dvar*100;

% Shrink the graph and position it
ah=gca;
ah.Position=ah.Position.*[1 3 1 0.5];

% Title the graph
title(sprintf('Variance reduction %i%% |\n Correlation coefficient %5.2f | Probability of chance %5.2f',...
              round(vardux),R,P))

% Print the graph
print -dpdf nasadata
