# Reproduce from the frozen NOAA file in the current directory.
# Base R only; no network calls or package installation.
stopifnot(unname(tools::md5sum('co2_mm_mlo.txt')) == '35591deeaf110bca96cc64a7ec337994')
d <- read.table('co2_mm_mlo.txt', comment.char='#', col.names=c('year','month','decimal','average','deseasonalized','days','sd','uncertainty'))
d <- subset(d, year <= 2025)
d$date <- as.Date(sprintf('%04d-%02d-01',d$year,d$month))
stopifnot(!anyDuplicated(d$date), identical(d$date,seq(min(d$date),max(d$date),by='month')))
d$flag_interpolated <- d$date >= as.Date('1974-05-01') & d$sd<0 & d$uncertainty<0
d$quality_ambiguous <- d$date >= as.Date('1974-05-01') & xor(d$sd<0,d$uncertainty<0)
d$used_co2 <- ifelse(d$average<0 | d$flag_interpolated,NA_real_,d$average)
write.csv(d,'monthly_clean.csv',row.names=FALSE,na='')
annual <- do.call(rbind,lapply(split(d,d$year),function(z) data.frame(year=z$year[1],n=sum(is.finite(z$used_co2)),mean_ppm=if(sum(is.finite(z$used_co2))==12) mean(z$used_co2) else NA_real_)))
annual$growth_ppm <- c(NA,diff(annual$mean_ppm))
write.csv(annual,'annual_summary.csv',row.names=FALSE,na='')
periods <- list('1960-1969'=1960:1969,'1970-1979'=1970:1979,'1980-1989'=1980:1989,'1990-1999'=1990:1999,'2000-2009'=2000:2009,'2010-2019'=2010:2019,'2020-2025'=2020:2025)
seasonal <- do.call(rbind,lapply(names(periods),function(p){
 z<-d[d$year %in% periods[[p]],]; z$centered_year<-z$decimal-mean(z$decimal)
 fit<-lm(used_co2~centered_year,data=z,na.action=na.exclude)
 z$departure<-residuals(fit)
 a<-aggregate(departure~month,z,mean)
 data.frame(period=p,month=a$month,departure_ppm=a$departure)
}))
growth <- do.call(rbind,lapply(names(periods),function(p){
 z<-annual[annual$year %in% periods[[p]],]
 s<-seasonal[seasonal$period==p,]
 data.frame(period=p,annual_mean_slope_ppm_per_year=unname(coef(lm(mean_ppm~year,z))[2]),mean_year_on_year_growth_ppm=mean(z$growth_ppm,na.rm=TRUE),seasonal_peak_to_trough_ppm=diff(range(s$departure_ppm)),peak_month=s$month[which.max(s$departure_ppm)],trough_month=s$month[which.min(s$departure_ppm)])
}))
write.csv(seasonal,'seasonal_profiles.csv',row.names=FALSE)
write.csv(growth,'period_summary.csv',row.names=FALSE)
png('co2_patterns.png',width=1800,height=1500,res=180)
par(mfrow=c(2,1),mar=c(4.3,4.5,3,1),las=1)
plot(d$date,d$used_co2,type='l',col='#777777',xlab='Year',ylab='Monthly CO2 (ppm)',main='Mauna Loa record: rising level and recurring seasonality')
lines(as.Date(paste0(annual$year,'-07-01')),annual$mean_ppm,col='#0072B2',lwd=2)
legend('topleft',c('Monthly record','Complete-year mean'),col=c('#777777','#0072B2'),lty=1,lwd=c(1,2),bty='n')
sel<-c('1960-1969','1990-1999','2010-2019','2020-2025')
cols<-c('#333333','#0072B2','#009E73','#D55E00')
plot(1:12,rep(0,12),type='n',ylim=range(seasonal$departure_ppm),xaxt='n',xlab='Calendar month',ylab='Detrended departure (ppm)',main='Seasonal profiles after within-period linear detrending')
abline(h=0,col='grey80')
axis(1,1:12,month.abb)
for(j in seq_along(sel)){z<-seasonal[seasonal$period==sel[j],];lines(z$month,z$departure_ppm,type='b',col=cols[j],lty=j,pch=c(16,17,15,18)[j],lwd=2)}
legend('bottomleft',sel,col=cols,lty=1:4,pch=c(16,17,15,18),bty='n',cex=.9)
dev.off()
cat('Annual endpoints:\n'); print(annual[annual$year %in% c(1959,2020,2025),])
print(growth)


y <- d$used_co2
origins <- which(d$date >= as.Date('2020-12-01') & d$date <= as.Date('2025-11-01'))
pred_list <- list(); diag_list <- list()
for(t in origins){
 training <- y[seq_len(t)]
 warnings_here <- character()
 fit <- tryCatch(withCallingHandlers(stats::arima(training,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),include.mean=FALSE,method='ML'),warning=function(w){warnings_here<<-c(warnings_here,conditionMessage(w));invokeRestart('muffleWarning')}),error=function(e)e)
 fit_ok <- !inherits(fit,'error') && fit$code==0
 fc <- if(fit_ok) tryCatch(predict(fit,n.ahead=12),error=function(e)NULL) else NULL
 i <- seq.int(max(13,t-59),t)
 pairs <- y[i]-y[i-12]
 drift <- if(sum(is.finite(pairs))>=12) mean(pairs,na.rm=TRUE) else NA_real_
 res <- if(fit_ok) as.numeric(residuals(fit)) else numeric()
 res <- res[seq_along(res)>13 & is.finite(res)]
 diag_list[[length(diag_list)+1]] <- data.frame(origin=as.character(d$date[t]),train_n=t,converged=fit_ok,code=if(inherits(fit,'error')) NA else fit$code,ma1=if(fit_ok)fit$coef[1] else NA,sma1=if(fit_ok)fit$coef[2] else NA,sigma=if(fit_ok)sqrt(fit$sigma2) else NA,ljung_box_24_p=if(length(res)>30)Box.test(res,lag=24,type='Ljung-Box',fitdf=2)$p.value else NA,acf1=if(length(res)>30)acf(res,plot=FALSE,lag.max=1)$acf[2] else NA,warnings=paste(warnings_here,collapse=';'),error=if(inherits(fit,'error'))conditionMessage(fit) else '')
 for(h in c(1,6,12)){
  j<-t+h
  if(j>nrow(d) || d$year[j]<2021 || d$year[j]>2025) next
  stopifnot(j-12<=t,t<j)
  p<-if(is.null(fc))NA_real_ else as.numeric(fc$pred[h])
  se<-if(is.null(fc))NA_real_ else as.numeric(fc$se[h])
  sn<-y[j-12]; sd<-sn+drift
  pred_list[[length(pred_list)+1]]<-data.frame(origin=as.character(d$date[t]),target=as.character(d$date[j]),horizon=h,train_n=t,actual=y[j],sarima=p,se=se,lo80=p-qnorm(.9)*se,hi80=p+qnorm(.9)*se,lo95=p-qnorm(.975)*se,hi95=p+qnorm(.975)*se,seasonal_naive=sn,seasonal_drift=sd,drift_annual=drift,drift_pairs=sum(is.finite(pairs)),target_observed=is.finite(y[j]),sarima_eligible=is.finite(p),naive_eligible=is.finite(sn),drift_eligible=is.finite(sd),common=all(is.finite(c(y[j],p,sn,sd))))
 }
}
predictions<-do.call(rbind,pred_list); diagnostics<-do.call(rbind,diag_list)
metrics<-do.call(rbind,lapply(c(1,6,12),function(h)do.call(rbind,lapply(c('sarima','seasonal_naive','seasonal_drift'),function(m){
 z<-predictions[predictions$horizon==h & predictions$common,]; e<-z[[m]]-z$actual
 data.frame(horizon=h,model=m,n=length(e),mae=mean(abs(e)),rmse=sqrt(mean(e^2)),bias=mean(e),coverage80=if(m=='sarima')mean(z$actual>=z$lo80 & z$actual<=z$hi80) else NA,coverage95=if(m=='sarima')mean(z$actual>=z$lo95 & z$actual<=z$hi95) else NA,width80=if(m=='sarima')mean(z$hi80-z$lo80) else NA,width95=if(m=='sarima')mean(z$hi95-z$lo95) else NA)
}))))
write.csv(predictions,'heldout_predictions.csv',row.names=FALSE,na='')
write.csv(metrics,'forecast_metrics.csv',row.names=FALSE,na='')
write.csv(diagnostics,'fit_diagnostics.csv',row.names=FALSE,na='')
png('forecast_evaluation.png',width=1800,height=1550,res=180)
par(mfrow=c(2,1),mar=c(4.3,4.5,3.1,1),las=1)
models<-c('sarima','seasonal_drift','seasonal_naive')
colors<-c('#0072B2','#009E73','#D55E00')
plot(c(1,6,12),rep(0,3),type='n',ylim=c(0,max(metrics$mae)*1.45),xlab='Forecast horizon (months)',ylab='Mean absolute error (ppm)',xaxt='n',main='Held-out targets in 2021-2025: expanding monthly refits')
axis(1,c(1,6,12))
for(k in seq_along(models)){z<-metrics[metrics$model==models[k],];lines(z$horizon,z$mae,type='b',lwd=2,pch=c(16,17,15)[k],lty=k,col=colors[k])}
legend('topleft',c('SARIMA','Seasonal naive + drift','Seasonal naive'),lty=1:3,pch=c(16,17,15),col=colors,bty='n')
z<-predictions[predictions$horizon==12,]; tt<-as.Date(z$target)
plot(tt,z$sarima-z$actual,type='n',ylim=range(z$lo95-z$actual,z$hi95-z$actual),xlab='Target month',ylab='Prediction - observation (ppm)',main='12-month SARIMA errors and conditional 95% intervals')
polygon(c(tt,rev(tt)),c(z$lo95-z$actual,rev(z$hi95-z$actual)),col='#C9E3F0',border=NA)
abline(h=0,lty=2,col='#444444');lines(tt,z$sarima-z$actual,col='#0072B2',lwd=2)
legend('bottomleft',c('Forecast error','95% interval relative to observed'),lty=c(1,NA),pch=c(NA,15),col=c('#0072B2','#C9E3F0'),bty='n',cex=.9)
dev.off()
print(metrics)
cat('Converged:',sum(diagnostics$converged),'/',nrow(diagnostics),'; warnings:',sum(nzchar(diagnostics$warnings)),'\n')
print(diagnostics[c(1,nrow(diagnostics)),])


# Independent boundary and arithmetic checks.
stopifnot(nrow(d)==814, max(d$date)==as.Date('2025-12-01'),
 all(as.Date(predictions$origin)<as.Date(predictions$target)),
 all(as.Date(predictions$target)>=as.Date('2021-01-01')),
 all(as.Date(predictions$target)<=as.Date('2025-12-01')),
 all(predictions$common),
 identical(as.integer(table(predictions$horizon)),c(60L,55L,49L)))
for(k in seq_len(nrow(predictions))){
 z<-predictions[k,]; t<-match(as.Date(z$origin),d$date); j<-match(as.Date(z$target),d$date)
 stopifnot(j-t==z$horizon, z$train_n==t, j-12<=t)
 stopifnot(abs(z$seasonal_naive-d$used_co2[j-12])<1e-10)
 ii<-seq(t-59,t)
 manual_drift<-mean(d$used_co2[ii]-d$used_co2[ii-12],na.rm=TRUE)
 stopifnot(abs(z$seasonal_drift-(d$used_co2[j-12]+manual_drift))<1e-10)
}
for(k in seq_len(nrow(metrics))){
 m<-metrics[k,]; z<-predictions[predictions$horizon==m$horizon & predictions$common,]
 stopifnot(abs(m$mae-sum(abs(z[[m$model]]-z$actual))/nrow(z))<1e-10)
}
checks<-data.frame(check=c('monthly grid and cutoff','forecast temporal boundaries','all baseline inputs at or before origin','baseline arithmetic','metric arithmetic','all 60 fits converged'),passed=c(TRUE,TRUE,TRUE,TRUE,TRUE,all(diagnostics$converged)))
write.csv(checks,'validation_checks.csv',row.names=FALSE)
capture.output(sessionInfo(),file='session_info.txt')
cat('Endpoint increase:',annual$mean_ppm[annual$year==2025]-annual$mean_ppm[annual$year==1959],'\n')
cat('Relative increase (%):',100*(annual$mean_ppm[annual$year==2025]/annual$mean_ppm[annual$year==1959]-1),'\n')
cat('Ljung Box p range:',range(diagnostics$ljung_box_24_p),'\n')
cat('Residual lag-1 ACF range:',range(diagnostics$acf1),'\n')
cat('Forecast-error lag-1 correlation (overlapping targets):\n')
for(h in c(1,6,12)){zz<-subset(predictions,horizon==h);cat(h,acf(zz$sarima-zz$actual,plot=FALSE,lag.max=1)$acf[2],'\n')}
print(checks)
