clear;
clc;

%==========================================================================
% MODEL DESCRIPTION
%
% X = st_ratio      
% Y = expend_pupil
% Z = score
% W = english proficiency
%
% Z = b0 + b1*X + b2*Y + b3*W + b4*X^2 + b5*W^2 + b6*X*Y
% 
% Fix W at some value, then
% Z = alpha + b1*X + b2*Y + b4*X^2 + b6*X*Y
%    where alpha = b0 + b3*W + b5*W^2
%==========================================================================

% parameter value setting
%------------------------------

b0 = 213.87;
b1 = 33.16;
b2 = 0.046;
b3 = -0.845;
b4 = -0.545;
b5 = 0.003;
b6 = -0.00227;

% fix W and define alpha 
%------------------------------
W = 50;
alpha = b0 + b3*W + b5*W^2;                 % alpha = b0 + b3*W + b5*W^2

% define mesh grid for X, Y 
%------------------------------

x = [0:0.25:30];                            % grid for X
y = [0:100:10000];                          % grid for Y

[X,Y] = meshgrid(x,y);                      % mesh grid for X, Y

% 3D - plot 
%------------------------------

Z = alpha + b1*X + b2*Y + b4*(X.*X) + b6*(X.*Y);    

figure(1)
surf(X,Y,Z);
xlabel('student-teacher ratio');
ylabel('expenditure per pupil');
zlabel('score');
