year = [1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022]
ΔT = [0.28, -0.19, 0.11, -0.02, 0.13, 0.16, 0.15, 0.33, 0.51, 0.14, 0.53, 0.3, 0.22, 0.31, 0.32, 0.56, 0.17, 0.36, 0.43, 0.46, 0.36, 0.27, 0.56, 0.25, 0.34, 0.6, 0.51, 0.34, 0.47, 0.71, 0.72, 0.61, 0.65, 0.5, 0.92, 0.27, 0.6, 0.73, 0.46, 0.44, 0.62, 0.69, 0.83, 1.12, 0.98, 0.75, 0.94, 1.14, 0.78, 0.89];

using PyPlot
plot(year, ΔT, "r.-")
xlabel("year")
ylabel("ΔT (°C) vs. 1901–2000 baseline")
title("Global average temperature change")

plot(year, ΔT, "r.-")
xlabel("year")
ylabel("ΔT (°C) vs. 1901–2000 baseline")
title("Global average temperature change")

for (i,j,c) in ((1,50,"k"), (26,34,"b"), (25,35,"g"))
    slope = (ΔT[j] - ΔT[i]) / (year[j] - year[i])
    plot(year[[i,j]], ΔT[[i,j]], "$(c)o")
    plot(year[[begin,end]], ΔT[i] .+ (year[[begin,end]] .- year[i]) .* slope, "$c-")
end

ylim(-0.2, 1.2)

b = ΔT
A = [ones(length(year)) year.-1973]

x̂ = A \ b

# Alternatively, the normal equations solution:
(A'A) \ A'b

using PyPlot
plot(year, ΔT, "r.-")
xlabel("year")
ylabel("ΔT (°C) vs. 1901–2000 baseline")
title("Global average temperature change: Linear fit")

plot(year, A * x̂, "k-")
legend(["data", "fit: slope $(round(x̂[2], sigdigits=2)) °C/year"])