How to calculate ARIMA(1,0,0)(1,0,1)12 prediction by hand

Let's simulate some data so we are discussing the same thing:

set.seed(1) # for reproducibility
foo <- ts(rnorm(120),frequency=12)
library(forecast)
model <- Arima(foo,order=c(1,0,0),seasonal=c(1,0,1))
summary(model)

This gives us coefficients as follows (among other information I cut out):

Coefficients:
          ar1     sar1    sma1    mean
      -0.0099  -0.1702  0.1455  0.1094
s.e.   0.0845   0.1488  0.1301  0.0778

The mean column gives the estimate of the intercept $c$, so $hat{c}=0.1094$. We also see that

$$ hat{phi}_1=-0.0099,quadhat{zeta}_1=-0.1702,quadhat{eta}_1=0.1455.$$

Now, don't confuse the standard errors of coefficients (the bottom row in that table) with the "errors" $e_t$ of the time series, which are often also called "innovations"! To calculate your forecast by hand, you will need the in-sample residuals, which you can get by residuals(model):

            Jan          Feb          Mar          Apr          May          Jun          Jul          Aug          Sep          Oct          Nov          Dec
1  -0.735553343  0.066955583 -0.943961915  1.476074374  0.234806723 -0.927359944  0.368693460  0.632512042  0.472470019 -0.409945615  1.397862585  0.294450877
2  -0.746077521 -2.329640213  0.969061751 -0.107573389 -0.121263921  0.810190805  0.729287954  0.507294672  0.826148959  0.670607748  0.006578943 -2.091766001
3   0.474295307 -0.218331743 -0.238847488 -1.591649628 -0.607182475  0.326703714  1.268850473 -0.189888368  0.294634012 -0.142091804 -1.493790693 -0.592080902
4  -0.494517692 -0.169220482  0.978319051  0.625789052 -0.281741566 -0.361444773  0.612571688  0.446748132 -0.789537517 -0.831412254  0.211147825  0.656072516
5  -0.229845305  0.764597566  0.322403594 -0.696611263  0.220069074 -1.246051099  1.321653144  1.896296747 -0.478259657 -1.177661559  0.460253138 -0.222673095
6   2.286638985 -0.106118277  0.582443469 -0.096566119 -0.847234754  0.041789417 -1.882604304  1.382027510  0.048955844  2.037863852  0.396060640 -0.824113175
7   0.550322552 -1.044477734 -1.359551587  0.169737483 -0.572862717 -0.107744203 -0.087922486 -0.672689316 -0.682331268 -0.196481153  1.074479552 -1.641277559
8   0.472271407  0.203579085  0.919979555 -0.400093431  0.246128733  0.156631224 -0.643675563  1.070879306  1.044609294  0.587086842  1.508505784  0.426411586
9  -1.370487392 -0.687123491 -1.311877423 -0.606575690 -0.727654495 -0.069989482 -1.038005727  0.068624068 -0.734687244  1.667274973  0.656806055  0.823737095
10  0.246994475  1.556966428 -0.766859421 -0.591604041  1.297917511 -0.749437420 -0.347055541 -0.508741319 -0.457409991 -0.354390605  0.391589413 -0.265387074

So, suppose we want to forecast for Jan 11. Your formula

$$ y_t- phi_1 y_{t-1} + zeta_1 phi_1 y_{t-13} - zeta_1 y_{t-12} = c + e_t - eta_1 e_{t-12}$$

turns into

$$ y_t= phi_1 y_{t-1} - zeta_1 phi_1 y_{t-13} + zeta_1 y_{t-12} + c + e_t - eta_1 e_{t-12}. $$

We replace the unknown parameters by their estimates as above. We take $y_{t-1}$, $y_{t-12}$ and $y_{t-13}$ from the series history, and take $e_{t-12}$ from the residuals(model) table - it's the entry for Jan 10, which is equal to $0.246994475$. Finally, we don't know $e_t$ yet, because it's our unknown new innovation, so we replace it by its expectation, which is zero. And there you are.

Note that I didn't calculate the actual predictions, because it's tedious, and because I am almost certain there is an error in the formula (but I still believe the description here is helpful). Please take a look at this earlier thread. Note in particular how Arima() with an intercept fits an ARIMA model to $Y_t-hat{c}$, so your formula will need to be adapted.