Fitting N datapoints in 3D on a straight line

I have N datapoints in 3d that lie on a line. The y-direction is fixed, so I want to fit x,z against y. Lets say we have 6 datapoints, that align with the y axis:

x=[0,0,0,0,0,0]
y=[1,2,3,4,5,6]
z=[0,0,0,0,0,0]

what I want to do:
I want to get the best set of fitting parameters, the gof and fitting error. So far with a least squarefit, I get a reduced chi2 of < 1, which means I might be overfitting (or misunderstanding something).
Questions:
1.) For example, for the above example I receive a reduced chi2 of 0- this seems false to me?
2.) Also, I am wondering if a least square fit is adequate for this as well- maybe someone can shed some insight on this? Would svd be a better choice for this?

import scipy.optimize
import numpy as np


#define a model (line)
def linear(params, y):
    a, b = params
    data = [a * y[i] + b for i in range(0, len(y))]
    return data

#define the residuals that need to me minimized
def fitting_cost(params, x, y, z):
    a_x, b_x, a_z, b_z = params
    x_pred = linear((a_x, b_x), y)
    z_pred = linear((a_z, b_z), y)
    res_x = [x_pred[i] - x[i] for i in range(0, 6)]
    res_z = [z_pred[i] - z[i] for i in range(0, 6)]
    return res_x + res_z


#do the fit and return parameters plus gof
def least_squares_fit(x, y, z):
    sp = [0,0,0,0]
    result = scipy.optimize.leastsq(fitting_cost, sp,
                                    args=(x, y, z),
                                    full_output=True)
    s_sq = (result[2]['fvec'] ** 2).sum() / (
            len(result[2]['fvec']) - len(result[0]))
    return result[0], s_sq