title: Running Generalized Least Squares in Python
body: <p>I am trying to do some regressions in Python using statsmodels.api, but my models all have problems with autocorrelation and heteroskedasticity. So I thought of trying out Generalized Least Squares (GLS). I am not very familiar with running this form of least squares, so stuck pretty close to the instructions on the below page:</p>
  
  <p><a href="https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.GLS.html" rel="nofollow noreferrer">https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.GLS.html</a></p>
  
  <p>with the below code:</p>
  
  <pre><code>X = b2_final_Train[explan_vars]
  Y = b2_final_Train[[ExognVar]]    
  model = sm.OLS(Y, X).fit()
  
  ##GLS Component    
  ols_resid = model.resid
  res_fit = sm.OLS(list(ols_resid[1:]), list(ols_resid[:-1])).fit()
  rho = res_fit.params
  
  from scipy.linalg import toeplitz
  order = toeplitz(np.arange(50014))
  sigma = rho**order
  
  gls_model = sm.GLS(Y, X, sigma=sigma)
  gls_results = gls_model.fit()
  print(gls_results.summary())
  </code></pre>
  
  <p>My problem as you can probably work out by looking at the code is sigma is a very big matrix, which overloads any computer I run it on being a 50014 x 50014 matrix. But as far as I am aware the GLS matrix is meant to be big enough for every error and that is how many observations I have in my data so that is how many errors I have. So is there something I am missing about running GLS which makes the problem computationally more manageable? </p>
  
  <p>Thanks </p>